Extended Geometry and Probability Model for GNSS + Constellation Performance Evaluation

: We proposed an extended geometry and probability model (EGAPM) to analyze the performance of various kinds of (Global Navigation Satellite System) GNSS + constellation design scenarios in terms of satellite visibility and dilution of precision (DOP) et al. on global and regional scales. Di ﬀ erent from conventional methods, requiring real or simulated satellite ephemerides, this new model only uses some basic parameters of one satellite constellation. Veriﬁed by the reference values derived from precise satellite ephemerides, the accuracy of visible satellite visibility estimation using EGAPM gets an accuracy better than 0.11 on average. Applying the EGAPM to evaluate the geometry distribution quality of the hybrid GNSS + constellation, where highly eccentric orbits (HEO), quasi-zenith orbit (QZO), inclined geosynchronous orbit (IGSO), geostationary earth orbit (GEO), medium earth orbit (MEO), and also low earth orbit (LEO) satellites included, we analyze the overall performance quantities of di ﬀ erent constellation conﬁgurations. Results show that QZO satellites perform slightly better in the Northern Hemisphere than IGSO satellites. HEO satellites can signiﬁcantly improve constellation geometry distribution quality in the high latitude regions. With 5 HEO satellites included in the third-generation BeiDou navigation satellite system (BDS-3), the average VDOP (vertical DOP) of the 30 ◦ N–90 ◦ N region can be decreased by 16.65%, meanwhile satellite visibility can be increased by 38.76%. What is more, the inclusion of the polar LEO constellation can signiﬁcantly improve GNSS service performance. When including with 288 LEO satellites, the overall DOPs (GDOP (geometric DOP), HDOP (horizontal DOP), PDOP (position DOP), TDOP (time DOP), and VDOP) are decreased by about 40%, and the satellite visibility can be increased by 183.99% relative to the Global Positioning System (GPS) constellation.


Introduction
Global Navigation Satellite System (GNSS) performance can be evaluated from several aspects, such as satellite signal quality [1,2], the precision of broadcast ephemerides [3][4][5], the multipath effect [6,7], convergence speed and precision of positioning [8,9]. Among all of these aspects, the satellite constellation, which defines the observation geometry distribution of navigation systems [10,11], is of fundamental importance, especially under the progress of constellation design. Satellite visibility and dilution of precision (DOP), which reflect the geometry distribution quality of constellations,

Methods
In the previous GAPM method [15,16], the observing sphere was divided into latitudes and longitudes grids by 1 • × 1 • , and the probability of a satellite appearing in one grid can be determined by its velocity. This probability is regarded as the satellite observing probability, which is the basis of the estimation of satellite visibility and DOP. For one specific station on the earth ellipsoid, the number of visible satellites is the sum of probability of the constellation in the observation area of the station. Furthermore, the satellite probability is regarded as weighting factors in the calculation of DOPs. Based on the DOPs, the performance and characterization of satellite visibility and geometry distribution for one GNSS constellation can be further evaluated. It is noted that the opportunities of a satellite's appearance in certain areas are zero where the latitude is higher than the satellite's orbit incline angle. Because the previous GAPM method is only applicable for the approximately circular orbit, we modify and extend the GAPM method in this paper to make it applicable for BDS-3 IGSO, QZO, and HEO satellites.

MEO and LEO Satellite Observing Probability
For GNSS MEO or LEO satellites in approximately circular orbits, previous GAPM is effective. The velocity of east-west direction V λ and north-south direction V ϕ for the position (ϕ,λ) can be written as [16]: where GM is the product of the gravitational constant and the earth's mass; R is the norm of a satellite position vector; i orb is the satellite's orbit incline angle to the equator; ϕ is geocentric latitude; λ is longitude.
As Equation (1) shows, the velocity of satellites appearing in one grid element is a function related to the latitude. So, for satellites appearing in the grid elements with the same geocentric latitude, the V ϕ term would be the same. Based on the hypothesis that GNSS satellites are distributed symmetrically in the east-west direction, the in and out satellites of the grid are balanced in this direction [15]. Therefore, the observing probability of a satellite depends on V ϕ only. The satellite observing probability can be modified as [15]: where k is the probability constant for a specific GNSS MEO or LEO satellite constellation. The superscript M/L stands for GNSS MEO or LEO satellites. As mentioned above, the spherical orbit surface can be divided into 1 • × 1 • grids centered in (ϕ i , λ j ): In addition, the probability constant k can be calculated as [16]: where n is the total number of GNSS MEO or LEO satellites.

GEO and IGSO Satellite Observing Probability
GEO satellites are relatively static to earth in earth-centered earth fixed (ECEF) system. Therefore, in theory, a GEO satellite will stay at the specific point all the time. For instance, the 3 BDS-3 GEO satellites are respectively located at 80 • E, 110.5 • E, and 140 • E [27]. So, in theory, the probability of a GEO satellite appearing at the designed position, such as 80 • E, is 1, while at other arbitrary positions is 0. Then GEO satellite observing probability can be written as [15]: For an IGSO satellite, their satellite ground track is a closed curve. For example, the 3 BDS-3 IGSO satellites, their ground tracks are coincident, while the longitude of the intersection point is Remote Sens. 2020, 12,2560 4 of 21 at 118 • E, with a phase difference of 120 • [28]. We can calculate a satellite position (X, Y, Z) in the ECEF system when the true anomaly f turns a specific angle and the set of these points in one earth self-rotation period is called B in this study. For an IGSO satellite, when f turns over from 0 • to 360 • , all coordinates in the ECEF system and corresponding velocity in one earth self-rotation period can be obtained. Considering both the accuracy and computation efficiency, this specific angle can be set at 0.25 • [15]. Then, refer to the idea of previous GAPM, the probability of IGSO satellites appearing at one point of B can be determined by the corresponding velocity. So, the modified probability definition, P IGSO (X,Y,Z) , of IGSO satellites to be observed at one point of B can be written as: where V IGSO (X,Y,Z) is the velocity at one point of B in the ECEF system; k IGSO is the probability constant for IGSO satellites which can be calculated by the following equation: where n IGSO is the total number of IGSO satellites with the same ground track. V IGSO (X,Y,Z) at one point of B can be calculated as follows. Firstly, a satellite position (X, Y, Z) in ECEF can be computed with the true anomaly f as the independent variable. Then, the velocity of X, Y, and Z components: V X , V Y , and V Z can be written as: where t is one specific epoch; X, Y, and Z are the ECEF coordinates of a satellite at the epoch t or t + ∆t; ∆t means the time elapsed from one point to another very adjacent point. Then, V (X,Y,Z) at one point is: We investigate the observing probability of the three BDS-3 IGSO satellites using the modified observing probability definition of IGSO satellites. From Figure 1, it can be seen that the observing probability is larger at the northern and southern ends of the figure-8-shape track which is closer to the actual situation than the corresponding results presented in [15].

HEO and QZO Satellite Observing Probability
As mentioned above, LEO, MEO, GEO, and IGSO satellites are all in approximately circle orbits with small eccentricity, and the magnitude of satellites' position vector and flight velocity vector are almost constant, which is the prerequisites for the earlier GAPM [15,16]. However, some special kinds of orbits have larger eccentricity, such as the QZO satellites of the Quasi-Zenith Satellite System (QZSS) [29]. The eccentricity of QZO satellites of about 0.075 is much larger than that of GPS satellites, which is normally less than 0.01. To investigate the performance of satellite orbits with large eccentricity, the traditional GAPM mentioned above should be generalized. For this purpose, the key issue is to calculate their space coordinates and corresponding velocity in the ECEF system. Therefore, some points need to be considered: (1) the magnitude of satellite position vector will change all the time; (2) the time that satellite fly from one point to another adjacent point is different because the satellite velocity and the magnitude of position vector are changing all the time; (3) the magnitude of satellite velocity in the ECEF system is changing also.

HEO and QZO Satellite Observing Probability
As mentioned above, LEO, MEO, GEO, and IGSO satellites are all in approximately circle orbits with small eccentricity, and the magnitude of satellites' position vector and flight velocity vector are almost constant, which is the prerequisites for the earlier GAPM [15,16]. However, some special kinds of orbits have larger eccentricity, such as the QZO satellites of the Quasi-Zenith Satellite System (QZSS) [29]. The eccentricity of QZO satellites of about 0.075 is much larger than that of GPS satellites, which is normally less than 0.01. To investigate the performance of satellite orbits with large eccentricity, the traditional GAPM mentioned above should be generalized. For this purpose, the key issue is to calculate their space coordinates and corresponding velocity in the ECEF system. Therefore, some points need to be considered: (1) the magnitude of satellite position vector will change all the time; (2) the time that satellite fly from one point to another adjacent point is different because the satellite velocity and the magnitude of position vector are changing all the time; (3) the magnitude of satellite velocity in the ECEF system is changing also.
For point (1), the magnitude of satellite position vector R can be calculated by the following equation [19]: where e is the orbit eccentricity; a is semi-major axis. From Equation (10), we know that the R is the function of the true anomaly f. For point (2), according to Kepler's second law, the radius vector sweeps over equal areas in equal time intervals [19], i.e., the area rate A  , is a constant: Therefore, the elapsed time for the satellite from one point to adjacent point can be calculated by the equation, A A /  , where A is the area that radius vector swept during the period. If the true anomaly changes by a very small amount and the area that the radius swept could be regarded as a For point (1), the magnitude of satellite position vector R can be calculated by the following equation [19]: where e is the orbit eccentricity; a is semi-major axis. From Equation (10), we know that the R is the function of the true anomaly f. For point (2), according to Kepler's second law, the radius vector sweeps over equal areas in equal time intervals [19], i.e., the area rate . A, is a constant: Therefore, the elapsed time for the satellite from one point to adjacent point can be calculated by the equation, A/ . A, where A is the area that radius vector swept during the period. If the true anomaly changes by a very small amount and the area that the radius swept could be regarded as a circular sector, then the area of this tiny circular sector, dA, can be calculated by the following equation: where df is the tiny change of true anomaly. The whole area could then be derived by the area integral within the interval [f, f + 0.25 • ]: However, the analytical expression of this integration is quite complex, and the result of numerical integration can be used. The interval [f, f + 0.25 • ] is divided into 10,000 pieces, and the above equation can be reformulated as: Remote Sens. 2020, 12, 2560 Thus, the elapsed time, ∆t, can be written as: When ∆t is determined, satellite velocity in the ECEF system could be derived. Satellite coordinates in the orbital plane coordinate system can be written as: where η is the axis points to the ascending node and ξ is the axis points to the direction on the argument of the latitude of 90 • . Satellite coordinates can be transformed from the orbital plane coordinate system into the ECEF system: where R 3 (−Ω i ) and R 1 (−i orb ) are rotation matrices as follows: where Ω i is the right ascension of ascending node which can be written as: where Ω 0 is the right ascending node at the initial epoch; ω earth is the angular velocity of the earth rotation. ∆t can be obtained by Equation (15). Thus, all ECEF cartesian coordinates in one earth rotation period can be calculated and the set of these coordinates is defined as Γ.
As for point (3), the magnitude of satellite velocity in the ECEF system can be calculated using the Equations (8) and (9).
Finally, the observing probability of a satellite at one point of Γ can be written as: where the superscript HEO/QZO stands for HEO or QZO satellites; k HEO/QZO is the corresponding probability constant for HEO or QZO satellites of large eccentricity orbits with the same ground track, which can be written as: is the velocity of HEO or QZO at one point of Γ in the ECEF system. Based on the mentioned above, observing the probability distribution of 3 QZO satellites can be seen from Figure 2. The eccentricity of QZO is set as 0.075, but the inclination and the center of the longitude of the ground track are the same as IGSO of BDS-3. The distribution shows that the probability is larger at the northern end of the track, which is very different from IGSO satellites, as a consequence of the velocity near apogee being lower. probability constant for HEO or QZO satellites of large eccentricity orbits with the same ground track, which can be written as: where / HEO QZO n is the total number of HEO or QZO satellites with the same ground track; is the velocity of HEO or QZO at one point of Γ in the ECEF system. Based on the mentioned above, observing the probability distribution of 3 QZO satellites can be seen from Figure 2. The eccentricity of QZO is set as 0.075, but the inclination and the center of the longitude of the ground track are the same as IGSO of BDS-3. The distribution shows that the probability is larger at the northern end of the track, which is very different from IGSO satellites, as a consequence of the velocity near apogee being lower.

Calculation of DOPs and Satellite Visibility
A satellite is visible when it is above a specific cutoff elevation angle (CEA). In this paper, to ensure higher accuracy, we apply the 5° CEA for the following DOP and satellite visibility calculations [15]. We calculated the probability for each satellite that the satellite elevation angle is higher than CEA for each grid element. Then the visible satellite number N can be written as:

Calculation of DOPs and Satellite Visibility
A satellite is visible when it is above a specific cutoff elevation angle (CEA). In this paper, to ensure higher accuracy, we apply the 5 • CEA for the following DOP and satellite visibility calculations [15]. We calculated the probability for each satellite that the satellite elevation angle is higher than CEA for each grid element. Then the visible satellite number N can be written as: The observation equation of a satellite can be written as [27]: where e x , e y , and e z are the unit vectors of station-to-satellite vector; v is the residual; cδt the receiver clock correction; l is the vector of observations; δX, δY, and δZ are the corrections to a station's approximate coordinates. In EGAPM, observation equations can be established for each grid or point higher than specific CEA and corresponding probability is regarded as its weight. Then, using the Remote Sens. 2020, 12, 2560 8 of 21 above information, for an arbitrary station on the earth ellipsoid, the coefficient matrix of the normal equations can be written as [16]: where D is the visible area above CEA. The inverse of the normal equation in Equation (25) gives the variance/covariance matrix Q: Then, the geometric dilution of precision, GDOP, by definition, is the square root of the trace of the variance/covariance matrix Q.
Cartesian coordinates can be translated into the station's topocentric up (U), north (N), and east (E) coordinate system using the following equations: where R 3 (−L) and R 2 (B) are rotation matrices which can be written as follow: Then, the normal equation can be written as: In addition, the covariance matrix Q UEN can be derived using the variance propagation rule as: Remote Sens. 2020, 12, 2560 So, by definition, the DOPs can be written as: Obtained from the above expression, the correlation, ρ Uclk , of station geodetic height between receiver clock difference follows the expression of: In addition, the ratio, NE, of the positioning precision in the longitude component to the latitude component can be written as: In summary, we can establish an observation equation based on all grids or points above CEA, and the probability of satellites in this grid is regarded as the weight of this observation equation. Then, the normal equation of the GNSS+ constellation can be given by: can be given by Equation (25).

Method Validation
We used 24 h averages on global scales of the values derived from real GPS precise ephemerides provided by GeoForschungsZentrum Potsdam (GFZ) and the precise ephemerides of the BDS-3+QZO simulation constellation, to validate the usefulness EGAPM and assess the method of satellite visibility and DOPs estimation. Constellation parameters of BDS-3+QZO and GPS are shown in Table 1. Table 1. Constellation parameters of third-generation BeiDou navigation satellite system (BDS-3)+ Quasi-Zenith orbit (QZO) and Global Positioning System (GPS) [28][29][30]. Estimations of the number of visible satellites and corresponding reference values at 5 • CEA are shown in Figure 3. Using EGAPM, the root mean squares (RMS) between the estimated and reference number of visible satellites are 0.10 for 32 GPS satellites and 0.05 for 30 BDS-3+QZO satellites. Generally, the discrepancies between estimated and reference values at high latitudes are larger than other locations.

BDS-3+QZO GPS
The underestimating rates for the GPS and BDS-3+QZO constellation are calculated (Table 2), to assess the DOPs underestimation with a 5 • CEA. The average underestimating rates are 7.95-11.07%. Estimations of the number of visible satellites and corresponding reference values at 5° CEA are shown in Figure 3. Using EGAPM, the root mean squares (RMS) between the estimated and reference number of visible satellites are 0.10 for 32 GPS satellites and 0.05 for 30 BDS-3+QZO satellites. Generally, the discrepancies between estimated and reference values at high latitudes are larger than other locations. The underestimating rates for the GPS and BDS-3+QZO constellation are calculated (Table 2), to assess the DOPs underestimation with a 5° CEA. The average underestimating rates are 7.95-11.07%.

Experimental Results and Analyses
In this section, the satellite visibility and geometry distribution quality of GNSS and GNSS+ constellations are investigated. Firstly, the performance of BDS-3 and GPS is compared. Then, 3 IGSO satellites of BDS-3 are replaced by 3 QZO satellites and positioning precision of IGSO and QZO satellites are compared. After that, 5 HEO satellites are combined with BDS-3 to establish the BDS-3+HEO constellation and the contribution of HEO is investigated. Finally, the performance of the GNSS+ constellation of GNSS+HEO+LEO is assessed.

BDS-3 VS GPS
Using the EGAPM the performance of the BDS-3 and GPS constellations is analyzed from the aspects of satellite visibility, DOPs, correlations level between station height and receiver clock difference, and the ratio of positioning precision of longitude and latitude components. The results are shown in Figure 4 and Table 3. GNSS+ constellation of GNSS+HEO+LEO is assessed.

BDS-3 VS GPS
Using the EGAPM the performance of the BDS-3 and GPS constellations is analyzed from the aspects of satellite visibility, DOPs, correlations level between station height and receiver clock difference, and the ratio of positioning precision of longitude and latitude components. The results are shown in Figure 4 and Table 3.    From Figure 4a,b, we can find that the satellite visibility of GPS is more globally balanced than BDS-3, and for BDS-3, more satellites can be observed in the Asia-Pacific region while fewer satellites in the Western Hemisphere. From Figure 4c to Figure 4l, the DOPs of BDS-3 are larger than GPS in the Western Hemisphere. Note that VDOP is larger in high latitude regions, especially in the polar regions for both BDS-3 and GPS because, with an increase in station latitude, the elevation of satellites observed becomes lower. In this case, introducing HEO satellites into the GNSS constellation may make more satellites be observed, which might improve the geometry distribution in the vertical direction. As shown in Figure 4m,n, for the GPS constellation in most regions, the positioning precision in the longitude component is worse than for the latitude component. However, in contrast, for the BDS-3 constellation in some low latitude regions and some polar regions, the precision of the longitude component is better. The reason may be that GPS satellites are all MEO satellites, mainly flying along east-west direction so positioning precision is better in this direction. However, in the case of BDS-3, 3 IGSO satellites mainly fly along the longitude direction. Besides, the GEO satellites mainly contribute to DOP value in the longitude direction mainly [2]. For the correlation level between station height and receiver clock difference, as shown in Figure 4o,p, BDS-3 in augmentation area (60 • S-60 • N and 50 • E-170 • E) is slightly higher than GPS in the same area. This may be because IGSO and GEO fly slowly and will persist over a user receiver for a longer time. Table 3 shows the average, minimal, and maximal visible satellite number and DOPs on global and regional scales. In the regional augmentation area, for BDS-3, 13.2 satellites can be seen on average, 2.7 and 2.2 more than the BDS-3 and GPS global average, respectively. Similarly, BDS-3 DOPs in the regional augmentation area is more favorable than BDS-3 and GPS on global scales. Because of fewer operational satellites of BDS-3 than GPS, on a global scale, visible satellites number and DOPs are slightly worse than GPS.

BDS-3+QZO Constellation
The QZSS with a service area covering the Asia-Pacific region is an augmentation system for GPS and Galileo [31]. The QZO satellites of QZSS are similar to the IGSO satellites of BDS-3. They all follow a figure-8-shape track on the spherical orbit surface in the ECEF system, but QZO satellite's ground track is not symmetric for the northern and southern hemisphere. This is because the eccentricity of QZO is 0.075 and therefore satellites with a higher elevation angle can be more easily observed in urban areas of Japan, which can improve the geometry distribution in the complex urban environment [32]. The specific parameters of QZO are listed in Table 4. Table 4. Quasi-Zenith orbit (QZO) parameters [29]. We replace 3 IGSO satellites of BDS-3 by 3 QZO satellites to investigate the differences between IGSO and QZO in terms of satellite visibility and GDOP. The inclination is set as 55 • and the center of longitude is modified as 118 • E which is the same as of BDS-3 IGSO satellites. The ground track of QZO satellites is shown in Figure 2. The values of the BDS-3+QZO constellation parameters are shown in Table 1.

Orbit Parameter Nominal Value
As revealed by Figure 5, after replacing 3 IGSO satellites, the performance in high latitude regions of the Northern Hemisphere is slightly improved but in some regions is slightly worse compared to BDS-3. This is because QZO is in an elliptical orbit and its apogee is in the Northern Hemisphere leading to QZO satellites spend about 2.3 h longer in the North Hemisphere.

BDS-3+HEO Constellation
Although the hybrid constellation is applied in BDS and QZSS to improve the visibility in some specific areas, the visibility in polar regions is still unfavorable. This is because that, for the current GNSS constellations, when the station is at higher latitudes, the elevation of satellites that can be observed is lower. If satellites that can cover high latitude regions are included in the GNSS+ constellation, the geometric distribution for such regions will be improved. The HEO satellites, for example, the Russian Molniya satellites, are the kind of satellites that can fulfill this function [19]. These satellites have synchronous 12 h orbits of 500-40,000 km altitude that are inclined at an angle of 63.4° to the equator which can not only ensure an exceptional coverage of the Northern Hemisphere but also minimizes the impact of orbital perturbations caused by the earth's oblateness. The HEO satellite's parameters applied in this study are listed in Table 5 and their ground track is continuously repeated as shown in Figure 6. The HEO in the 3D graphic view is displayed in Figure  7. Table 5. Highly eccentric orbitsHEO) parameters.

Orbit Parameter
Nominal Value

BDS-3+HEO Constellation
Although the hybrid constellation is applied in BDS and QZSS to improve the visibility in some specific areas, the visibility in polar regions is still unfavorable. This is because that, for the current GNSS constellations, when the station is at higher latitudes, the elevation of satellites that can be observed is lower. If satellites that can cover high latitude regions are included in the GNSS+ constellation, the geometric distribution for such regions will be improved. The HEO satellites, for example, the Russian Molniya satellites, are the kind of satellites that can fulfill this function [19]. These satellites have synchronous 12 h orbits of 500-40,000 km altitude that are inclined at an angle of 63.4 • to the equator which can not only ensure an exceptional coverage of the Northern Hemisphere but also minimizes the impact of orbital perturbations caused by the earth's oblateness. The HEO satellite's parameters applied in this study are listed in Table 5 and their ground track is continuously repeated as shown in Figure 6. The HEO in the 3D graphic view is displayed in Figure 7.   Five HEO satellites and BDS-3 satellites are combined to establish the BDS-3+HEO constellation. The primary constellation parameters are listed in Table 6. We evaluate the performance of the BDS-3+HEO constellation on a global scale, as illustrated in Figure 8.  Five HEO satellites and BDS-3 satellites are combined to establish the BDS-3+HEO constellation. The primary constellation parameters are listed in Table 6. We evaluate the performance of the BDS-3+HEO constellation on a global scale, as illustrated in Figure 8. Five HEO satellites and BDS-3 satellites are combined to establish the BDS-3+HEO constellation. The primary constellation parameters are listed in Table 6. We evaluate the performance of the BDS-3+HEO constellation on a global scale, as illustrated in Figure 8. The detailed improvement rates of introducing HEO satellites relative to BDS-3 at 30° N-90° N region are summarized in Table 7. As can be seen from Figure 8 and Table 7, the DOPs and satellite visibility are modified favorably after using 5 HEO satellites. In the 30° N-90° N region, the satellite visibility is improved on average by 38.76%. In addition, all kinds of DOPs are improved. By introducing five HEO satellites to BDS-3, on average, values of the GDOP, PDOP, HDOP, VDOP, and TDOP are reduced by 15.70%, 16.65%, 10.35%, 16.65%, and 11.52% respectively. However, the correlation level between station height and receiver clock difference is still higher. As shown in Figure 8h, the correlation level in low latitude and the polar region is slightly higher. What is more, as revealed in Figure 8g and Figure 9, after applying HEO satellites, in some regions such as north polar and low latitude regions, positioning precision in longitude is better than in the latitude component. This may be since IGSO and HEO satellites mainly fly along the north-south direction, which would contribute to the DOP value in the longitude component mainly. The detailed improvement rates of introducing HEO satellites relative to BDS-3 at 30 • N-90 • N region are summarized in Table 7. As can be seen from Figure 8 and Table 7, the DOPs and satellite visibility are modified favorably after using 5 HEO satellites. In the 30 • N-90 • N region, the satellite visibility is improved on average by 38.76%. In addition, all kinds of DOPs are improved. By introducing five HEO satellites to BDS-3, on average, values of the GDOP, PDOP, HDOP, VDOP, and TDOP are reduced by 15.70%, 16.65%, 10.35%, 16.65%, and 11.52% respectively. However, the correlation level between station height and receiver clock difference is still higher. As shown in Figure 8h, the correlation level in low latitude and the polar region is slightly higher. What is more, as revealed in Figures 8g  and 9, after applying HEO satellites, in some regions such as north polar and low latitude regions, positioning precision in longitude is better than in the latitude component. This may be since IGSO and HEO satellites mainly fly along the north-south direction, which would contribute to the DOP value in the longitude component mainly.

BDS-3+LEO+HEO Constellation
The augmentation of GNSS with LEO satellites is one of the most interesting trends in the future. We consider the addition of 288 LEO satellites in polar orbits which is mentioned in Li et al. [25] to form the BDS-3+LEO+HEO constellation. The primary constellation parameters are listed in Table 8. These LEO orbits are distributed in 12 planes and orbit altitudes of 1000 km as shown in Figure 10 in the 3D graphic view. Using the EGAPM, the performance of the BDS-3+LEO+HEO constellation is analyzed. The results are plotted in Figure 11.

BDS-3+LEO+HEO Constellation
The augmentation of GNSS with LEO satellites is one of the most interesting trends in the future. We consider the addition of 288 LEO satellites in polar orbits which is mentioned in Li et al. [25] to form the BDS-3+LEO+HEO constellation. The primary constellation parameters are listed in Table 8. These LEO orbits are distributed in 12 planes and orbit altitudes of 1000 km as shown in Figure 10 in the 3D graphic view. Using the EGAPM, the performance of the BDS-3+LEO+HEO constellation is analyzed. The results are plotted in Figure 11.

BDS-3+LEO+HEO Constellation
The augmentation of GNSS with LEO satellites is one of the most interesting trends in the future. We consider the addition of 288 LEO satellites in polar orbits which is mentioned in Li et al. [25] to form the BDS-3+LEO+HEO constellation. The primary constellation parameters are listed in Table 8. These LEO orbits are distributed in 12 planes and orbit altitudes of 1000 km as shown in Figure 10 in the 3D graphic view. Using the EGAPM, the performance of the BDS-3+LEO+HEO constellation is analyzed. The results are plotted in Figure 11.   As can be seen in Figure 11, after introducing the LEO satellites, the DOPs have decreased significantly, which means that the geometric distribution has been improved. Since all LEO satellites are in polar orbits, a lot more satellites can be observed in high latitude regions than in the middle and low latitude areas (Figure 11a). Although hundreds of LEO satellites in polar orbits are introduced, the correlation between station height and receiver clock difference is still higher ( Figure  11h). As for the positioning precision ratio between the longitude and latitude, illustrated in Figure  11g, Figure 9, and Figure 8g. We can find that LEO satellites mainly contribute to longitude component performance. Table 9 reveals the average, minimal, and maximal improving rates of visible satellites and DOPs on a global scale. More satellites, with improving rate of 198% and 184%, respectively, for BDS-3 and GPS, can be observed. Besides, DOP values are about 40% lower than BDS-3 and GPS, which may improve the standard point positioning (SPP) accuracy and accelerate the precise point positioning (PPP) convergence speed.  As can be seen in Figure 11, after introducing the LEO satellites, the DOPs have decreased significantly, which means that the geometric distribution has been improved. Since all LEO satellites are in polar orbits, a lot more satellites can be observed in high latitude regions than in the middle and low latitude areas (Figure 11a). Although hundreds of LEO satellites in polar orbits are introduced, the correlation between station height and receiver clock difference is still higher (Figure 11h). As for the positioning precision ratio between the longitude and latitude, illustrated in Figure 11g, Figure 9, and Figure 8g. We can find that LEO satellites mainly contribute to longitude component performance. Table 9 reveals the average, minimal, and maximal improving rates of visible satellites and DOPs on a global scale. More satellites, with improving rate of 198% and 184%, respectively, for BDS-3 and GPS, can be observed. Besides, DOP values are about 40% lower than BDS-3 and GPS, which may improve the standard point positioning (SPP) accuracy and accelerate the precise point positioning (PPP) convergence speed.

Discussion
This article aims to provide a convenient way of assessing the global distribution of DOP and satellite visibility for the arbitrary navigation satellite constellation design. To do so, we replace the conventional point-by-point simulation of the constellation by statistical considerations of the probability of finding in certain viewing directions. We extend an earlier work of Wang et al. [15] and Wang et al. [16] to include non-circular orbits in analysis and present results for various sets of mixed-orbit constellations. Unfortunately, a disadvantage is also existing that EGAPM is not suitable to assess the ambiguity dilution of precision (ADOP), which can predict when ambiguity resolution can be expected to be successful [31,[33][34][35]. Applying EGAPM, if a large number of grids or points can be viewed regarded as satellites and their ambiguities are estimated, the computation burden will be too heavy. What is worse, the calculation results of ADOP are not coincident with corresponding results calculated by GFZ precise ephemerides. However, in our next research, we will commit to finding a more ingenious way to solve this problem.

Conclusions
Hybrid constellations will be one of the development trends in the future; hence, it is necessary to develop a simple strategy to evaluate the performance of GNSS+ constellations before these satellites are launched. We extend the previous GAPM, which is only suitable for circular orbits. Furthermore, most calculations of satellite visibility and DOPs are based on real or complex simulated ephemerides. However, our strategy estimates satellite visibility and DOPs based on a few orbital parameters of one satellite constellation.
The features of the BDS-3 and GPS constellation were investigated. Satellite visibility of GPS is more balanced on a global scale than BDS-3. In contrast, for BDS-3, more satellites can be observed in the Asia-Pacific region and fewer satellites in the Western Hemisphere. VDOP is larger, especially in the polar regions for both BDS-3 and GPS. Furthermore, it is interesting to find that for the BDS-3 constellation, the precision in the longitude component is better in some low latitude and some polar regions, which are in contrast to GPS in these areas. Besides, the correlation level between station height and receiver clock difference is higher for both BDS-3 and GPS.
For QZO satellites, their orbit's eccentricity is 0.075, and apogee is located above the Northern Hemisphere. Thus, these satellites stay a longer time above the Northern Hemisphere. Experiments show that when the BDS IGSO satellite is replaced by the QZO satellite, satellite visibility and GDOP are slightly better in the Northern Hemisphere but worse in the Southern Hemisphere relative to the BDS-3 constellation.
HEO satellites are introduced to improve satellite visibility in high latitude regions. The performance of the BDS-3+HEO constellation is assessed, and results show that the geometric distribution in the middle and high latitude regions is significantly improved. It is worth noting that VDOP is decreased by 16.65%, and satellite visibility improved by 38.76% on average in the 30 • N-90 • N region.
The performance of LEO satellites is also investigated, and the BDS-3+HEO+LEO constellation is analyzed. By introducing hundreds of LEO satellites, DOPs decreased greatly, 40% lower or so relative to the BDS-3 and GPS constellation on a global scale. In the meantime, it is interesting to find that the ratio of positioning precision between the longitude and the latitude component is strongly related to the direction of satellite travel. As LEO satellites are all at polar orbits, mainly flying along north-south direction and experiment results reveal that positioning accuracy is better in the longitude component in most regions.
The results and conclusion of this paper present a preliminary exploration for the design of the GNSS+ constellation and thus may benefit future GNSS development.