Integrated Precise Orbit Determination of Multi-GNSS and Large LEO Constellations

: Global navigation satellite system (GNSS) orbits are traditionally determined by observation data of ground stations, which usually need even global distribution to ensure adequate observation geometry strength. However, good tracking geometry cannot be achieved for all GNSS satellites due to many factors, such as limited ground stations and special stationary characteristics for the geostationary Earth orbit (GEO) satellites in the BeiDou constellation. Fortunately, the onboard observations from low earth orbiters (LEO) can be an important supplement to overcome the weakness in tracking geometry. In this contribution, we perform large LEO constellation-augmented multi-GNSS precise orbit determination (POD) based on simulated GNSS observations. Six LEO constellations with di ﬀ erent satellites numbers, orbit types, and altitudes, as well as global and regional ground networks, are designed to assess the inﬂuence of di ﬀ erent tracking conﬁgurations on the integrated POD. Then, onboard and ground-based GNSS observations are simulated, without regard to the observations between LEO satellites and ground stations. The results show that compared with ground-based POD, a remarkable accuracy improvement of over 70% can be observed for all GNSS satellites when the entire LEO constellation is introduced. Particularly, BDS GEO satellites can obtain centimeter-level orbits, with the largest accuracy improvement being 98%. Compared with the 60-LEO and 66-LEO schemes, the 96-LEO scheme yields an improvement in orbit accuracy of about 1 cm for GEO satellites and 1 mm for other satellites because of the increase of LEO satellites, but leading to a steep rise in the computational time. In terms of the orbital types, the sun-synchronous-orbiting constellation can yield a better tracking geometry for GNSS satellites and a stronger augmentation than the polar-orbiting constellation. As for the LEO altitude, there are almost no large-orbit accuracy di ﬀ erences among the 600, 1000, and 1400 km schemes except for BDS GEO satellites. Furthermore, the GNSS orbit is found to have less dependence on ground stations when incorporating a large number of LEO. The orbit accuracy of the integrated POD with 8 global stations is almost comparable to the result of integrated POD with a denser global network of 65 stations. In addition, we also present an analysis concerning the integrated POD with a partial LEO constellation. The result demonstrates that introducing part of a LEO constellation can be an e ﬀ ective way to balance the conﬂict between the orbit accuracy and computational e ﬃ ciency. plays a dominant role in non-gravitational for LEO. All the POD computations are performed in the single-thread mode on HP Appollo2000 machines with 16-core CPU and


Introduction
The precise orbit and clock products of global navigation satellite system (GNSS) are of great importance for GNSS applications, such as precise point positioning (PPP), GNSS meteorology etc. using simulated data. Their initial results demonstrated the significant accuracy improvements of GPS and BDS orbits as a result of the LEO constellation.
The aforementioned works evidence the contribution of onboard GNSS observations to the GNSS POD, and shed light on the performance of the integrated POD. However, these studies mainly focus on the integrated POD with a few LEO satellites, since no onboard GNSS observations from large LEO constellations are available at present. Additionally, the impact of different LEO configurations, such as LEO number, orbital altitude, and orbit type, on integrated POD has not yet been well studied. Motivated by the aforementioned studies, we propose a LEO constellation-augmented multi-GNSS precise orbit determination method and assess the LEO-augmented POD performance with different configurations of ground stations and LEO constellations. Within this paper, we mainly focus on the geometric factors that limit the accuracy of GNSS orbit determination. Six LEO constellations with different orbital altitudes, satellite numbers, and orbit types are designed to evaluate the influence of LEO constellations on the GNSS POD. Several station networks with different station number and distributions are employed to discuss the dependency of the integrated POD on ground stations. The paper is structured as follows. Section 2 starts by introducing GNSS and LEO constellations we used, and gives a detailed description of the sophisticated simulation processing methods used for GNSS observations. Then, the LEO constellation-augmented four-system precise orbit determination algorithm and strategy are introduced in Section 3. Section 4 analyzes LEO-augmented GNSS POD results under different tracking conditions in detail. A discussion of our result is given in Section 5. Finally, conclusions are provided in Section 6.

Constellation Design
With the full operation of the Galileo and BDS constellations, there will be more than 120 GNSS satellites in the sky by 2020. As of August 2019, Galileo had reached its nominal constellation, with only back-up satellites to be launched in the future, while the BDS constellation, which is still under construction, has 19 BDS-3 satellites providing navigation services. Therefore, in order to perform the integrated POD with the whole constellation of the four systems, the orbits of the GPS, GLONASS, Galileo, and BDS satellites have to be simulated based on the nominal constellation configurations [22][23][24][25] of the four systems. The nominal GPS constellation consists of 24 medium earth orbit (MEO) satellites unevenly distributed in six planes. MEO satellites are also employed by GLONASS and Galileo using a 24/3/1 Walker constellation to provide global coverage. In addition, six additional MEO satellites serve as in-orbit spare satellites for the Galileo constellation. Different from the other three constellations, BDS-3 is made up of 24 MEO satellites, 3 geostationary earth orbit (GEO) satellites, and 3 inclined geosynchronous orbit (IGSO) satellites. The three GEO satellites are located at 80 • E, 110.5 • E, and 140 • E, while the IGSO satellites operate in an orbit with an inclination of 55 • ; the right ascension of the ascending nodes (RAAN) is 118 • E. Table 1 presents detailed information on the GNSS constellations. To comprehensively evaluate the impact of LEO configuration on GNSS POD, six LEO constellations with different configurations are designed and simulated. Table 2 lists the LEO constellation parameters. Nearly-polar and sun-synchronous orbits are selected in our study to investigate the influence of the LEO orbit type, which are two typical orbit types for many LEO missions. The nearly-polar orbit we choose has the same orbit inclination, i.e., 84.6 degrees, as the Iridium satellite constellation. Considering the fact that computation time rises sharply for the integrated POD due to huge amounts of observations to be processed and unknown parameters to be recovered, only the constellations with satellites number of 60, 66, 96 are adopted. The Chinese Hongyan project, which has been announced by the China Aerospace Science and Technology Corporation (CASC) for the purpose of global communication and LEO-augmented positioning, proposes a constellation made up of 60 LEO satellites as its backbone, while the 66-LEO constellation is employed by the Iridium satellite constellation for global coverage. In order to assess the influence of LEO altitude on the LEO-augmented GNSS POD, the 60-LEO constellation is simulated with altitudes of 600 km, 1000 km, and 1400 km. Figure 1 presents a sketch of the designed LEO constellations.

Observations Simulation Configuration
At present, no onboard multi-GNSS observations from large LEO constellations are available, because all the LEO constellations providing augmented navigation services are still under construction. Hence, in order to achieve an augmented four-system POD, all kinds of data should be simulated. Since we mainly focus on the contribution of onboard multi-GNSS observations to the GNSS POD, only ground and onboard multi-GNSS measurements are taken into consideration in the

Observations Simulation Configuration
At present, no onboard multi-GNSS observations from large LEO constellations are available, because all the LEO constellations providing augmented navigation services are still under construction.
Hence, in order to achieve an augmented four-system POD, all kinds of data should be simulated. Since we mainly focus on the contribution of onboard multi-GNSS observations to the GNSS POD, only ground and onboard multi-GNSS measurements are taken into consideration in the data simulation.
In this study, we adopt the same simulation method described in Li et al. [9]. To make the simulated data as close as possible to the real observations, all errors related to the satellite, receiver, and signal path, as well as observation noise, should be considered and calculated in the process of data simulation. The undifferenced code and carrier phase observations for ground stations and LEO satellites can be expressed as follows: where s, g, leo, and j represent GNSS satellites, ground stations, LEO satellites and signal frequency respectively. P and L denote pseudorange and carrier phase observations respectively. ρ s g is the distance between the mass center of satellites and ground receiver, while ρ s leo is the distance between the mass center of both GNSS satellites and LEO satellites. c is the speed of light in a vacuum. δt g , δt leo , and δt s represent the clock offsets for the ground receiver, onboard receiver, and satellite respectively. I s g,j and I s leo,j are the ionospheric delays at frequency j. T s g is the tropospheric delay of ground station. λ j is the wavelength at the frequency j. N s g,j and N s leo,j refers to the integer ambiguities for ground and onboard receiver respectively. (b g, j , b leo,j , b s j ) and (B g,j , B leo, j , B s j ) represent the code hardware delay and carrier phase hardware delay of the ground receiver, onboard receiver, and satellites respectively. (ε s g,j , ε s leo,j ) and (ω s g,j , ω s leo,j ) are the combination of multipath and noise for code and carrier phase observations. The main task of the GNSS observation simulations is to accurately compute the components on the right side of Equation (1) using the existing models, and to guarantee that the simulated GNSS observations reflect a real-world environment as much as possible. Before the observation simulation, the orbits of both GNSS and LEO satellites are firstly simulated. Then, we performed a standard multi-GNSS PPP using the real GNSS measurements of each ground station to provide a weekly solution of ground station position, receiver clock offset, the wet components of tropospheric delay, inter-system bias (ISB), and inter-frequency bias (IFB) for the subsequent observation simulation. The satellite-to-receiver distance is calculated based on the position for station at the signal receiving time and mass center position for satellites at the signal transmitting time. The previous PPP processing provides a weekly solution of the ground station position, while the GNSS satellites position can be obtained from the simulated orbits. The phase center offset (PCO) and phase center variation (PCV) values of satellites and ground stations should be considered using the values from "igs14.atx" [26], though they are not presented in Equation (1). The multi-GNSS receiver clock offsets are simulated using the PPP-derived receiver clock offsets, as well as both ISB and IFB. The multi-GNSS precise clock products from GeoForschungsZentrum (GFZ) are used to provide the values of satellite clock offsets. Given the ongoing development of BDS and Galileo constellations, the clock offsets of the unoperated satellites are replaced with the existing one, e.g., the clock offset of BDS C27 satellite is considered the same as that of C13 satellite. It's noteworthy that both the satellite precise clock products and PPP-derived receiver clock offsets are generated based on an ionosphere-free (IF) combination, and they usually absorb the IF combination of code hardware delays for the satellite and receiver, respectively, Remote Sens. 2019, 11, 2514 6 of 22 due to the strong coupling between the clock offsets and the code hardware delay. Hence, Equation (1) can be rewritten as with The combination of code hardware delays in the code observation equation can be calculated using known differential code biases (DCB) from Center for Orbit Determination in Europe (CODE) [27,28]. The ionospheric delay along the direction of the signal propagation at a specific frequency can be modeled using the international GNSS service (IGS) global ionosphere maps (GIM). The dry component of slant tropospheric delay is computed using the Saastamoinen model [29], as well as the global mapping function (GMF) [30], while the wet component is derived as a parameter of the PPP solution. We set the carrier phase ambiguities of each continuous arc as integer values and the phase delay was assumed to be a small floating constant. The residuals of code and carrier phase observations of the PPP process are employed to calculate multipath errors and noise corrections. Moreover, the relativistic correction, phase wind-up correction, and tidal displacements are also taken into consideration in the simulation.
In the LEO onboard observations simulation, the calculation method of the components associated with GNSS satellites is the same as that for the ground measurements simulation. Similar to GNSS satellites, the position of the LEO mass center can be acquired from the simulated LEO orbits. The PCO and PCV values of onboard antennas are set to zero. The receiver clock offsets of ground stations are used with the ISB and IFB values derived from GFZ multi-GNSS bias products to simulate the multi-GNSS clock offsets for each LEO onboard receiver. Different from the ground stations, the signals between the GNSS and LEO satellites only pass through the topside part of the ionosphere. Hence, the ionospheric delay for onboard observations is generated according to its contribution to total electron content computed from GIM [31]. In terms of observation noise, the standard deviations of the random noise for code and phase observations are set to 1 m and 5 mm respectively. Details of the employed models and simulation standards are described in Li et al. [9].

Integrated POD Method
Based on the simulated GNSS observations, the LEO-augmented four-system POD, also known as the integrated POD, can be performed. The ionosphere-free (IF) combination is adopted to eliminate the ionospheric delay in the dual-frequency code and phase measurements. The linearized observation equation for IF combinations can be written as follows: O leo,0 = x leo,0 , y leo,0 , z leo,0 , v leo,x , v leo,y , v leo,z , p leo,1 , p leo,2 , · · · , p leo,n where u s g and u s leo denote the unit vector of the direction from ground receiver to satellite and from LEO to satellite respectively. ϕ(t, t 0 ) s and ϕ(t, t 0 ) leo refer to the state-transition matrix for satellite and LEO, which transfer the satellites state from reference epoch t 0 to epoch t. O s 0 and O leo,0 are the initial state vector for GNSS and LEO satellites respectively, which consist of the position and velocity of the satellite at the initial epoch and dynamics parameters. For GNSS satellites, the estimated dynamics parameters mainly refer to the solar radiation pressure (SRP) parameters, while for LEO satellites, these dynamics parameters are usually made up of SRP parameters, atmosphere drag parameters, and empirical accelerations. M s g is the mapping function for tropospheric delay, and Z g denotes the zenith delay of the tropospheric wet component. N s g,IF and N s leo,IF represent the float ambiguity for ground station and LEO satellites respectively. The other symbols have similar meanings as those described in Equation (2), except for the IF linear combination.
In the case of the multi-GNSS integrated POD, the estimated parameters can be divided into satellite-dependent parameters, ground station-dependent parameters, and LEO-dependent parameters, which can be expressed as follows: where superscript G, R, E, and C represent GPS, GLONASS, Galileo, and BDS satellites respectively.

Processing Strategy
We select about 120 globally-distributed stations from the Multi-GNSS Experiment (MGEX) [1], and 22 global stations from the International GNSS Monitoring and Assessment System (iGMAS) [32], to simulate one week of ground GNSS observations from day of year (DOY) 001-007, 2018. Figure 2 indicates the distribution of ground stations we used. The onboard GNSS measurements are also simulated based on the LEO constellations described in Section 2.1. The simulated observation interval is set to 30 s for both ground and onboard observations. Table 3 lists the processing strategy of integrated multi-GNSS POD in detail. Due to the large amount of observations to be processed, we selected an arc length of 24 h and a processing interval of 300 s. The cut-off elevation angle is set to 7 • and 1 • for ground stations and LEO satellites respectively. In terms of force model, GNSS and LEO satellites suffer from different perturbative forces, since they move at different orbital altitudes, especially in the aspect of non-gravitational forces. For GNSS satellites, the solar radiation pressure serves as the primary source of non-gravitational forces due to the thin atmosphere at the GNSS altitude; the atmosphere drag is neglected in the GNSS POD processing. Different from GNSS satellites, as the trajectories of LEO satellites undergo perturbation from both solar radiation pressure and atmosphere drag, the atmosphere drag plays a dominant role in the non-gravitational forces for LEO. All the POD computations are performed in the single-thread mode on HP Appollo2000 machines with 16-core CPU and 96-GB memory.

Processing Strategy
We select about 120 globally-distributed stations from the Multi-GNSS Experiment (MGEX) [1], and 22 global stations from the International GNSS Monitoring and Assessment System (iGMAS) [32], to simulate one week of ground GNSS observations from day of year (DOY) 001-007, 2018. Figure 2 indicates the distribution of ground stations we used. The onboard GNSS measurements are also simulated based on the LEO constellations described in Section 2.1. The simulated observation interval is set to 30 s for both ground and onboard observations.  Table 3 lists the processing strategy of integrated multi-GNSS POD in detail. Due to the large amount of observations to be processed, we selected an arc length of 24 h and a processing interval of 300 s. The cut-off elevation angle is set to 7° and 1° for ground stations and LEO satellites respectively. In terms of force model, GNSS and LEO satellites suffer from different perturbative forces, since they move at different orbital altitudes, especially in the aspect of non-gravitational forces. For GNSS satellites, the solar radiation pressure serves as the primary source of nongravitational forces due to the thin atmosphere at the GNSS altitude; the atmosphere drag is neglected in the GNSS POD processing. Different from GNSS satellites, as the trajectories of LEO satellites undergo perturbation from both solar radiation pressure and atmosphere drag, the atmosphere drag plays a dominant role in the non-gravitational forces for LEO. All the POD computations are performed in the single-thread mode on HP Appollo2000 machines with 16-core CPU and 96-GB memory. Table 3. Detailed processing strategy for the integrated multi-GNSS POD.

Result Analysis and Discussion
In this section, the performance of the multi-GNSS integrated POD with different LEO constellations and different station network schemes is analyzed. Several impact factors including LEO numbers, orbital types, and altitude, as well as the numbers and distribution of ground stations, are discussed in detail. The previously simulated orbits are regarded as the true values, and the differences between our estimated orbits and true orbits provide a method by which to assess the accuracy of the integrated POD results.

Integrated POD with Different Numbers of LEO
We firstly performed the integrated POD using global network from MGEX. The ground observations from the 120-station MGEX network are processed as a reference. In order to make a fair comparison, only a 65-station network with a good global distribution from the MGEX we selected is used to perform the integrated POD (shown as red dots with a black circle edge in Figure 2). The onboard observations from 60, 66, and 96 polar-orbiting LEO satellites with altitudes of 1000 km (see Figure 1a,b,d) are processed to investigate the influence of LEO numbers on the integrated POD. Figures 3-6 present the average 3D Root Mean Square (RMS) values of orbit differences between integrated POD solution and true orbits for all satellites from GPS, GLONASS, Galileo, and BDS constellations respectively. The result shows that by using only the ground multi-GNSS observations, the majority of GNSS satellites can achieve a 3D orbit accuracy of better than 4 cm, which is comparable to the accuracy of the 24-h orbit recovered from the real ground tracking measurements. This demonstrates that the error models and simulation strategy are close to the real situation. Table 4 computes the average 3D RMS of orbit differences for four systems. The worst orbit quality can be found in BDS GEO, and the average 3D RMS values of orbit differences even exceeds 2 m. This is reasonable, because the stationary and regional coverage feature of GEO satellites lead to poor tracking geometry, greatly hampering the orbit determination for GEO satellites. Galileo, and BDS constellations respectively. The result shows that by using only the ground multi-GNSS observations, the majority of GNSS satellites can achieve a 3D orbit accuracy of better than 4 cm, which is comparable to the accuracy of the 24-h orbit recovered from the real ground tracking measurements. This demonstrates that the error models and simulation strategy are close to the real situation. Table 4 computes the average 3D RMS of orbit differences for four systems. The worst orbit quality can be found in BDS GEO, and the average 3D RMS values of orbit differences even exceeds 2 m. This is reasonable, because the stationary and regional coverage feature of GEO satellites lead to poor tracking geometry, greatly hampering the orbit determination for GEO satellites.       As shown in Figures 3-6, all satellites providing global coverage can achieve a sub-centimeter level orbit with the support of LEO onboard observations, while the orbit accuracy of regional coverage satellites (BDS GEO and IGSO) is at the centimeter level. It can be seen that LEO satellites can make more contributions to GNSS orbits than the same numbers of ground stations, since they can not only provide high-quality observations, but also improve the geometry diversity. The orbit of GNSS satellites obtains a remarkable accuracy improvement after introducing the LEO observations into the POD processing. The accuracy improvement with respect to the ground-based POD result can reach over 70% for all of 60 LEO-, 66 LEO-, and 96 LEO-augmented POD schemes (as shown in Table 4). This indicates that the onboard data from the large LEO constellations can be an important supplement to the GNSS POD, even under the condition of a global network. We find that the BDS GEO satellites present the largest accuracy improvement when the integrated POD is implemented; the improvement percentage even reaches up to about 98%. This inspiring result evidences the contribution of LEO to improve the tracking geometry of GEO satellites.
Furthermore, it can be seen that smaller orbit differences can be achieved when more LEO satellites are introduced into the integrated POD. The 96-LEO solution achieves a slightly better orbit accuracy than the other two LEO-augmented POD solutions. However, compared with the slight accuracy improvement of the 96-LEO solution, introducing additional 30 LEO satellites into the processing of the integrated POD leads to a huge burden on the computational efficiency because of the larger number of unknown parameters to be solved. The average computation time of integrated POD for the 96-LEO scheme is about 38 h, which is almost three times more than that of 60-LEO and 66-LEO schemes (about 13 h). As shown in Figures 3-6, all satellites providing global coverage can achieve a sub-centimeter level orbit with the support of LEO onboard observations, while the orbit accuracy of regional coverage satellites (BDS GEO and IGSO) is at the centimeter level. It can be seen that LEO satellites can make more contributions to GNSS orbits than the same numbers of ground stations, since they can not only provide high-quality observations, but also improve the geometry diversity. The orbit of GNSS satellites obtains a remarkable accuracy improvement after introducing the LEO observations into the POD processing. The accuracy improvement with respect to the ground-based POD result can reach over 70% for all of 60 LEO-, 66 LEO-, and 96 LEO-augmented POD schemes (as shown in Table 4). This indicates that the onboard data from the large LEO constellations can be an important supplement to the GNSS POD, even under the condition of a global network. We find that the BDS GEO satellites present the largest accuracy improvement when the integrated POD is implemented; the improvement percentage even reaches up to about 98%. This inspiring result evidences the contribution of LEO to improve the tracking geometry of GEO satellites.
Furthermore, it can be seen that smaller orbit differences can be achieved when more LEO satellites are introduced into the integrated POD. The 96-LEO solution achieves a slightly better orbit accuracy than the other two LEO-augmented POD solutions. However, compared with the slight accuracy improvement of the 96-LEO solution, introducing additional 30 LEO satellites into the processing of the integrated POD leads to a huge burden on the computational efficiency because of the larger number of unknown parameters to be solved. The average computation time of integrated POD for the 96-LEO scheme is about 38 h, which is almost three times more than that of 60-LEO and 66-LEO schemes (about 13 h). To assess the impact of LEO orbital type, two typical LEO orbits, nearly-polar and sun-synchronous orbits, are adopted in this study. Both of LEO constellations consist of 60 LEO satellites flying at an altitude of 1000 km (see Figure 1a,b). The ground observations from a global network with 65 stations are used. The results are shown in Figure 7. It can be seen that the sun-synchronous-orbiting constellation presents a slightly stronger enhancement to the GNSS orbits than the polar-orbiting constellation. With the sun-synchronous-orbiting data of 60 satellites, the orbit of GNSS satellites can achieve an accuracy similar to that of the 96-polar-orbiting scheme, and the orbit accuracy of BDS GEO for the sun-synchronous scheme is even better than that for the 96-polar-orbiting scheme.  To assess the impact of LEO orbital type, two typical LEO orbits, nearly-polar and sunsynchronous orbits, are adopted in this study. Both of LEO constellations consist of 60 LEO satellites flying at an altitude of 1000 km (see Figure 1 (a) and (b)). The ground observations from a global network with 65 stations are used. The results are shown in Figure 7. It can be seen that the sunsynchronous-orbiting constellation presents a slightly stronger enhancement to the GNSS orbits than the polar-orbiting constellation. With the sun-synchronous-orbiting data of 60 satellites, the orbit of GNSS satellites can achieve an accuracy similar to that of the 96-polar-orbiting scheme, and the orbit accuracy of BDS GEO for the sun-synchronous scheme is even better than that for the 96-polarorbiting scheme. As reported by Li et al. [9], the polar-orbiting constellation outperforms the sun-synchronousorbiting constellation in terms of the convergence time of PPP. However, the situation is reversed As reported by Li et al. [9], the polar-orbiting constellation outperforms the sun-synchronous-orbiting constellation in terms of the convergence time of PPP. However, the situation is reversed when the LEO constellation is employed to enhance the GNSS precise orbit determination. In order to determine the reason for this, we computed the actual number of LEO satellites used for the GNSS POD and calculated an orbit dilution of precision (ODOP) value for every GNSS satellite using a similar algorithm like the position dilution of precision (PDOP) for ground stations. Assuming one GNSS satellite can be tracked by n ground stations simultaneously, the ODOP of this GNSS satellite can be calculated as follows: where (x n , y n , z n ) and (x s , y s , z s ) are the positions of ground station and GNSS satellite respectively. → r n = (x n − x s , y n − y s , z n − z s ) is the distance vector with the direction from GNSS satellite to ground station. The ODOP value describes the geometric strength of observations for GNSS satellites. The better the geometric strength, the smaller the ODOP value. Figure 8 shows the used LEO numbers and ODOP values for C02 (GEO) and G13 (MEO) respectively. The two 60-LEO schemes contribute an almost comparable number of LEO satellites to the orbit determination for BDS GEO C02. The average number of the used sun-synchronous-orbiting LEO satellite is 21.5, which is slightly larger than polar-orbiting LEO number, i.e., 20.6. In terms of ODOP, the ODOP values of the sun-synchronous-orbiting constellation present a stronger fluctuation with a larger amplitude than those of polar-orbiting constellation. For BDS GEO satellites, the changeless tracking geometry is the key limitation of their orbit determination. The main contribution of LEO satellites is to bring a fast variation to the tracking geometry for BDS GEO. The faster fluctuation of ODOP indicates that the sun-synchronous-orbiting constellation can bring more variations to the GEO tracking geometry, which can yield more pronounced orbit accuracy improvements for GEO satellites. In the case of satellite G13, the differences between the ODOP variation of the two 60-LEO solutions is small, which leads to a relatively small discrepancy in the corresponding orbit accuracy, i.e., 1 mm.   Table 5 lists the result of the integrated POD with LEO satellites at different altitudes. The 60 polar-orbiting LEO satellites at altitudes of 600, 1000, and 1400 km are considered (see Figure 1(a)). It can be seen that there are almost no orbit accuracy differences among the solutions. The integrated POD with higher LEO satellites presents a slightly higher orbit quality than that with lower LEO satellites, except for the BDS GEO satellites. For BDS GEO, the best orbit accuracy is achieved by the 600 km solution. Given the motionless characteristic of GEO satellites, the LEO satellites orbiting at  Table 5 lists the result of the integrated POD with LEO satellites at different altitudes. The 60 polar-orbiting LEO satellites at altitudes of 600, 1000, and 1400 km are considered (see Figure 1a). It can be seen that there are almost no orbit accuracy differences among the solutions. The integrated POD with higher LEO satellites presents a slightly higher orbit quality than that with lower LEO satellites, except for the BDS GEO satellites. For BDS GEO, the best orbit accuracy is achieved by the 600 km solution. Given the motionless characteristic of GEO satellites, the LEO satellites orbiting at an altitude of 600 km can contribute more improvements to the tracking geometry for BDS GEO due to their faster motion, which may be the reason for the better performance of the 600-km solution in GEO orbit accuracy. In addition, it should be noted that the LEO satellites at low altitude usually suffer from complex SRP modelling. Also, the trajectory of low-altitude LEO is perturbed by a stronger drag force, which will increase the difficulty of drag force modeling and lead to accelerated orbital decay for LEO. The difficulty of air drag and SRP modelling may negatively affect the result of the integrated POD. Fortunately, with the proper modeling of these forces, such as estimating the drag parameter with a shorter interval, low-altitude LEO satellites can obtain a high precision orbit result, which has been evidenced in many low-altitude LEO missions, such as Swarm [38].

Integrated POD with Different Ground Network
In order to evaluate the impact of station quantity and distribution on the LEO-augmented POD, fewer global stations are selected to perform the integrated POD with different LEO constellations. We firstly adopted about 22 global stations from iGMAS (see Figure 2). Typically, the three LEO schemes, i.e., 60 polar-orbiting LEO, 60 sun-synchronous-orbiting LEO, and 96 polar-orbiting LEO constellations at the altitude of 1000 km, are chosen. Figure 9 illustrates the integrated POD results for GPS, GLONASS, Galileo, and BDS using iGMAS ground observations. For comparison, we also plot the GNSS orbit results determined from the 65 MGEX stations and 60 sun-synchronous-orbiting LEO satellites. The corresponding statistical values are provided in Table 6. With regard to the numbers provided in Table 4, a clear degradation in the orbit accuracy can be recognized for all GNSS satellites compared with the POD result using MGEX network when only ground iGMAS observations are processed. This is reasonable, because fewer ground stations are employed to recover GNSS satellites orbit. Similar to the MGEX solution, the significant orbit accuracy improvement for all satellites can be recognized after the incorporation of LEO onboard observations into the GNSS orbit determination. Certainly, BDS GEO satellites exhibit the largest orbit accuracy improvement. The different LEO constellations present a similar performance to that in the MGEX solution. As shown in Table 6, the sun-synchronous-orbiting constellation yields slightly better performance in LEO-augmented POD compared to the two polar-orbiting constellations.
fewer ground stations are employed to recover GNSS satellites orbit. Similar to the MGEX solution, the significant orbit accuracy improvement for all satellites can be recognized after the incorporation of LEO onboard observations into the GNSS orbit determination. Certainly, BDS GEO satellites exhibit the largest orbit accuracy improvement. The different LEO constellations present a similar performance to that in the MGEX solution. As shown in Table 6, the sun-synchronous-orbiting constellation yields slightly better performance in LEO-augmented POD compared to the two polarorbiting constellations. Figure 9. Average 3D orbit differences RMS of GNSS satellites for the integrated POD using the iGMAS network. P and S denote the nearly-polar and sun-synchronous orbits, respectively. Figure 9. Average 3D orbit differences RMS of GNSS satellites for the integrated POD using the iGMAS network. P and S denote the nearly-polar and sun-synchronous orbits, respectively. Meanwhile, it should be noted that although a sparse iGMAS network is employed to perform the integrated POD, the orbit accuracy of all three iGMAS integrated POD schemes is comparable to that of the corresponding MGEX solutions. The performance of the iGMAS and MGEX solutions mainly differs in their estimated BDS GEO orbits. Compared with the MGEX solution, the integrated POD with the iGMAS network provides a relatively worse orbit for the BDS GEO. The result demonstrates that the dependence of GNSS orbit determination on the ground stations can be largely reduced when introducing a large number of LEO satellites into the POD processing. Serving as moving stations, LEO satellites can not only provide a large amount of onboard observations, but also bring evident improvements to the tracking geometry for GNSS satellites.
The above results indicate that the onboard observations play an important role in reducing the contribution of ground observations to GNSS orbit determination when introducing a large number of LEO satellites. To further investigate the dependency of the integrated POD on the ground stations, we design three integrated POD schemes with only 8 or 4 ground stations involved. The distribution of the small number of stations we selected is shown in Figure 10. The sun-synchronous-orbiting constellation with 60 LEO satellites is adopted due to its best performance in the previous study. It should be noted that the earth rotation parameters are fixed when we perform the integrated POD with a few stations.
contribution of ground observations to GNSS orbit determination when introducing a large number of LEO satellites. To further investigate the dependency of the integrated POD on the ground stations, we design three integrated POD schemes with only 8 or 4 ground stations involved. The distribution of the small number of stations we selected is shown in Figure 10. The sun-synchronous-orbiting constellation with 60 LEO satellites is adopted due to its best performance in the previous study. It should be noted that the earth rotation parameters are fixed when we perform the integrated POD with a few stations.  As shown in Figure 11, the estimated orbit of the integrated POD with 8 regional stations can achieve a centimeter-level accuracy, which is better than the POD result only using ground MGEX observations (see Table 4). This indicates that with the assistance of 60 LEO satellites, only using observations from regional stations can already obtain a relatively high-quality satellite orbit. A significant reduction in orbit difference can be recognized for all the GNSS satellites when the 8 regional stations are replaced by 8 global stations. We can find that the orbit accuracy of the integrated POD with 8 global stations is almost in line with the accuracy of the corresponding MGEX and iGMAS integrated POD schemes. The accuracy improvement can be attributed to the better performance of the globally distributed stations in anchoring the whole constellations compared to the regional stations. Slight orbit degradation for GNSS satellites is observed when the number of global stations is reduced to 4. The result shows that the POD of GNSS satellites has less dependence on the number of ground stations and is more sensitive to the distribution of ground stations after the inclusion of the large LEO constellation. As shown in Figure 11, the estimated orbit of the integrated POD with 8 regional stations can achieve a centimeter-level accuracy, which is better than the POD result only using ground MGEX observations (see Table 4). This indicates that with the assistance of 60 LEO satellites, only using observations from regional stations can already obtain a relatively high-quality satellite orbit. A significant reduction in orbit difference can be recognized for all the GNSS satellites when the 8 regional stations are replaced by 8 global stations. We can find that the orbit accuracy of the integrated POD with 8 global stations is almost in line with the accuracy of the corresponding MGEX and iGMAS integrated POD schemes. The accuracy improvement can be attributed to the better performance of the globally distributed stations in anchoring the whole constellations compared to the regional stations. Slight orbit degradation for GNSS satellites is observed when the number of global stations is reduced to 4. The result shows that the POD of GNSS satellites has less dependence on the number of ground stations and is more sensitive to the distribution of ground stations after the inclusion of the large LEO constellation.

Integrated POD with Partial LEO Constellation
The previous results demonstrate the strong enhancement of the entire LEO constellation to the GNSS orbit estimation. However, it can be seen that the integrated POD suffers from a long computation time when a large number of LEO satellites are introduced. The computation time of

Integrated POD with Partial LEO Constellation
The previous results demonstrate the strong enhancement of the entire LEO constellation to the GNSS orbit estimation. However, it can be seen that the integrated POD suffers from a long computation time when a large number of LEO satellites are introduced. The computation time of the 60-LEO scheme is about 13 h, while for 96-LEO scheme, it reaches about 38 h. Only moderate accuracy differences are obtained for the 60-, 66-and 96-LEO constellations. In order to balance the conflict between orbit accuracy and computational efficiency, it is worthwhile to investigate the integrated POD with a partial LEO constellation, and to explore how many LEO satellites can be used as an auxiliary to improve the orbit accuracy of GNSS satellites. In this section, we process the observation from 65 MGEX stations and a partial LEO constellation. The 60 polar-orbiting LEO constellation and 60 sun-synchronous-orbiting LEO constellation at an altitude of 1000 km are chosen. Figure 12 shows the average 3D RMS of GNSS orbit differences for the integrated POD as a function of the number of LEOs. We find that with only adding 10 LEO satellites, BDS GEO satellite orbits can already achieve an accuracy of better than 10 cm, while the orbit differences of other satellites are below 2 cm. The corresponding computation time is less than 2 h (shown in Figure 13). This result indicates that the introduction of a small number of LEO satellites can not only obtain a relatively high accuracy GNSS orbit, but that it also limits the computation time, which is very beneficial to the GNSS application with high timeliness requirements. For both, the orbit differences of GNSS satellites decrease monotonically as the number of introduced LEO satellites increases. However, with more LEO satellites included, the improvement of orbit accuracy gradually becomes smaller. A reduction of about 0.5-4 cm can be observed in the orbit differences when the LEO increases from 10 to 20, whereas the additional 10 LEO satellites can only contribute to an accuracy improvement of less than 0.2 cm as the LEO increases from 40 to 50. On the other hand, the increase of LEO satellites introduces a large number of parameters to be estimated, resulting in a steep rise in computational time (shown in Figure 13). For example, the computation time of the integrated POD greatly rise from 9.46 h to 13.42 h when the LEO number increases from 50 to 60. But the addition of 10 LEO satellites only yields an orbit accuracy improvement of about 1 mm for the majority of GNSS satellites. satellites are below 2 cm. The corresponding computation time is less than 2 h (shown in Figure 13). This result indicates that the introduction of a small number of LEO satellites can not only obtain a relatively high accuracy GNSS orbit, but that it also limits the computation time, which is very beneficial to the GNSS application with high timeliness requirements. For both, the orbit differences of GNSS satellites decrease monotonically as the number of introduced LEO satellites increases. However, with more LEO satellites included, the improvement of orbit accuracy gradually becomes smaller. A reduction of about 0.5-4 cm can be observed in the orbit differences when the LEO increases from 10 to 20, whereas the additional 10 LEO satellites can only contribute to an accuracy improvement of less than 0.2 cm as the LEO increases from 40 to 50. On the other hand, the increase of LEO satellites introduces a large number of parameters to be estimated, resulting in a steep rise in computational time (shown in Figure 13). For example, the computation time of the integrated POD greatly rise from 9.46 h to 13.42 h when the LEO number increases from 50 to 60. But the addition of 10 LEO satellites only yields an orbit accuracy improvement of about 1 mm for the majority of GNSS satellites. In addition, it can be seen that the sun-synchronous LEO satellites exhibit a stronger enhancement to the GNSS orbits than those in nearly-polar orbit, no matter how many LEO satellites are introduced. The accuracy differences between the sun-synchronous LEO scheme and nearly-polar LEO scheme is very small for the GPS, GLONASS, Galileo, BDS IGSO, and MEO satellites. This is because the relative motion between these satellites and ground stations already yields good geometric diversity, so that the small contribution differences in tracking geometry between sun- In addition, it can be seen that the sun-synchronous LEO satellites exhibit a stronger enhancement to the GNSS orbits than those in nearly-polar orbit, no matter how many LEO satellites are introduced. The accuracy differences between the sun-synchronous LEO scheme and nearly-polar LEO scheme is very small for the GPS, GLONASS, Galileo, BDS IGSO, and MEO satellites. This is because the relative motion between these satellites and ground stations already yields good geometric diversity, so that the small contribution differences in tracking geometry between sun-synchronous LEO and nearly-polar LEO satellites have little impact on the orbit determination of these satellites. However, differences of more than 1 cm are visible in the BDS GEO orbits between the two solutions. Evidently, the combination of global stations and sun-synchronous LEO satellites result in better quality GEO orbits. This is because LEO satellites in sun-synchronous orbit can give rise to faster variations of observation geometry for BDS GEO satellites (see Figure 8).

Discussion
Serving as moving stations, the LEO satellites can significantly strengthen the tracking geometric diversity of GNSS satellites and improve the accuracy of the GNSS precise orbit determination, particularly in the case of regional network or sparse global network. Although many studies have reported the contribution of LEO onboard observations to the GNSS POD, few studies have focused on the performance of the integrated POD with a large LEO constellation. Meanwhile, the impact factors of the integrated POD, such as LEO number, orbital height, and orbit type, is barely analyzed. In this study, we presented the GNSS orbit solution derived in an integrated processing of the ground network and the large LEO constellation using simulated data. Six LEO constellations are designed to study the performance of the integrated POD with different tracking configurations. The results demonstrate that the orbit accuracy of GNSS satellites can be dramatically improved when a large number of LEO satellites is introduced; the more LEO satellites, the better the orbit accuracy. Notably, BDS GEO orbits present an accuracy of a few centimeters with the recruitment of LEO satellites. Moreover, the LEO satellite in sun-synchronous orbit can contribute more to GNSS orbits than nearly polar-orbiting satellites, since they give rise to a stronger observation geometry. Furthermore, the orbit altitude of LEO satellites is found to have no evident impact on the GNSS satellites except BDS GEO satellites. In addition, it can be seen that the integrated POD with the entire LEO constellation suffers from a long computational time because of the abundant GNSS observations to be processed and the large number of parameters to be estimated. The long computation time of the integrated POD cannot meet the timeliness requirement of the real-time precise positioning. This weakness in terms of computational efficiency can be overcome by just introducing a certain number of LEO satellites.
Although our study presents an optimistic integrated POD result due to the simulated data, the impact of introducing LEO constellations to improve GNSS orbit determination is clearly demonstrated. Indeed, not only tracking geometry, but many other factors, such as attitude model, solar radiation pressure and drag models, and antenna calibrations, constitute key limitations in realworld LEO and GNSS POD. Fortunately, with the efforts of many scientists, the force and observation models have been rigorously refined in recent years, which can reduce the impact of model errors on the GNSS POD as much as possible. Using the current models, GNSS satellites (except BDS GEO) can

Discussion
Serving as moving stations, the LEO satellites can significantly strengthen the tracking geometric diversity of GNSS satellites and improve the accuracy of the GNSS precise orbit determination, particularly in the case of regional network or sparse global network. Although many studies have reported the contribution of LEO onboard observations to the GNSS POD, few studies have focused on the performance of the integrated POD with a large LEO constellation. Meanwhile, the impact factors of the integrated POD, such as LEO number, orbital height, and orbit type, is barely analyzed. In this study, we presented the GNSS orbit solution derived in an integrated processing of the ground network and the large LEO constellation using simulated data. Six LEO constellations are designed to study the performance of the integrated POD with different tracking configurations. The results demonstrate that the orbit accuracy of GNSS satellites can be dramatically improved when a large number of LEO satellites is introduced; the more LEO satellites, the better the orbit accuracy. Notably, BDS GEO orbits present an accuracy of a few centimeters with the recruitment of LEO satellites. Moreover, the LEO satellite in sun-synchronous orbit can contribute more to GNSS orbits than nearly polar-orbiting satellites, since they give rise to a stronger observation geometry. Furthermore, the orbit altitude of LEO satellites is found to have no evident impact on the GNSS satellites except BDS GEO satellites. In addition, it can be seen that the integrated POD with the entire LEO constellation suffers from a long computational time because of the abundant GNSS observations to be processed and the large number of parameters to be estimated. The long computation time of the integrated POD cannot meet the timeliness requirement of the real-time precise positioning. This weakness in terms of computational efficiency can be overcome by just introducing a certain number of LEO satellites.
Although our study presents an optimistic integrated POD result due to the simulated data, the impact of introducing LEO constellations to improve GNSS orbit determination is clearly demonstrated. Indeed, not only tracking geometry, but many other factors, such as attitude model, solar radiation pressure and drag models, and antenna calibrations, constitute key limitations in real-world LEO and GNSS POD. Fortunately, with the efforts of many scientists, the force and observation models have been rigorously refined in recent years, which can reduce the impact of model errors on the GNSS POD as much as possible. Using the current models, GNSS satellites (except BDS GEO) can achieve centimeter-level orbits [1], though more efforts are needed to further refine the force and observation models. Meanwhile, the force model errors do not represent a weakness in the observability and tracking geometry, which is a separate issue from the tracking geometry that needs to be further investigated. In addition, except for the factors mentioned above, there are still many issues which need to be taken into consideration when implementing the integrated POD using real GNSS observations. For example, the incorporation of a large number of real multi-GNSS onboard observations may introduce unexpected unmodeled errors into the orbit recovery process. Meanwhile, different from ground receivers, the performance of space-borne receivers can be more easily affected by the space environment due to the high altitude, i.e., 600 km~1000 km, where LEO satellites are flying. As discussed in Xiong et al. [39], large equatorial plasma irregularities can degrade the tracking capability of the Swarm onboard receivers, leading to severe signal loss at low latitudes. This loss of signal for space-borne receivers undoubtedly has negative effects on the integrated POD. Moreover, except for GNSS data, many other types of satellite tracking data are not considered in our study. Once the construction of the navigation augmentation LEO constellation is complete, the new ranging observations from LEO satellites to ground stations can be employed. These new types of satellite tracking data are expected to remove the potential systemic bias in the estimated orbits and further improve the orbit accuracy of both GNSS and LEO satellites.

Conclusions
This paper is devoted to investigating the integrated precise orbit determination with large LEO constellations. We performed LEO-augmented multi-GNSS precise orbit determination in this study. Based on simulated GNSS observations, several integrated POD schemes are implemented to investigate the performance of integrated POD under different tracking conditions. The potential influence factors of the integrated POD, including LEO satellites number, altitude, and orbit type as well as ground station number and distribution, are analyzed and discussed in detail.
The result shows that joint processing ground observations from the global network and the large LEO constellation onboard observations can significantly improve the orbit accuracy for GNSS satellites, especially for BDS GEO satellites. The accuracy improvement percentage with respect to the ground-based POD results can reach over 70% for all the integrated POD schemes with 60-, 66-, 96-LEO satellites. The largest orbit accuracy improvement of over 98% can be recognized for BDS GEO, since the fast motion of LEO satellites brings tremendous variations to the tracking geometry of GEO satellites. By incorporating a large number of LEO satellites, BDS GEO satellites can obtain a centimeter-level orbit. Compared with the 60-and 66-LEO schemes, a slightly better orbit quality is observed in the 96 LEO scheme due to the introduction of more LEO satellites. However, the increase of the involved LEO satellites results in a sharp rise in the computational time of the integrated POD because of more unknown parameters to be solved. We also present the impact of the LEO orbit type on the integrated POD. With the same number of LEO satellites, the sun-synchronous-orbiting constellation presents a stronger enhancement to GNSS orbits than the polar-orbiting constellation. The improvement can be attributed to the more rapid variation of GNSS tracking geometry provided by a sun-synchronous-orbiting constellation. In terms of the LEO altitude, the orbit altitude of LEO satellites is found to have little influence on the enhancement to GNSS orbits except for with BDS GEO satellites. Benefiting from the faster motion of lower LEO satellites, the 600 km scheme achieves a better orbit accuracy for BDS GEO than the 1000 km and 1400 km schemes.
In order to assess the impact of station number and distribution, global and regional networks with different numbers of stations are employed. With the inclusion of large LEO constellation, the integrated POD appears to be more sensitive to the distribution of ground stations, rather than the station number. This is because in the integrated POD, the dependence of GNSS orbit determination on ground stations is largely reduced by the onboard GNSS observations when a large LEO constellation is considered. Based on 8 regional stations, the orbit accuracy of the integrated POD with a sun-synchronous-orbiting constellation is at the centimeter level. However, a subcentimeter accuracy can be recognized for the orbits of GNSS satellites when 8 globally distributed stations are adopted, which is comparable to the result of the integrated POD with a denser global network of 65 stations.
Although GNSS orbit estimations can greatly benefit from the incorporation of a full LEO constellation, the integrated POD suffers from long computation time due to the addition of a large number of LEO satellites. To balance the conflict between orbit accuracy and computational efficiency, the integrated POD with a partial LEO constellation is investigated. The orbit accuracy of GNSS satellites improves gradually as the LEO satellites increase. The more LEO satellites, the better the orbit accuracy. However, with an increasing number of LEOs, the accuracy improvement of GNSS orbits becomes smaller. Meanwhile, the increase of LEO satellites results in a steep rise in computational time. Considering both orbit accuracy and computational efficiency, we prefer 40 LEO satellites in a sun-synchronous orbit.
With the rapid development of the LEO constellations, we expect an improvement of GNSS POD when jointly processing GNSS observation data from reference sites and LEOs. Continued efforts for force modeling, algorithm optimization, and scheme design are required to achieve better performance of the integrated POD.