1. INTRODUCTION
One of the major differences of the Galileo satellite navigation system compared to the existing Global Positioning System (GPS) is the embedded integrity function, which enables the Galileo system to provide navigation for safety critical applications such as civil aviation. In the current Galileo baseline, the integrity aspects concerning the Signal-In-Space Error (SISE), which is the satellite-to-user error due to satellite clock and orbit errors, will be achieved by two parameters: Signal-In-Space Accuracy (SISA) and Integrity Flag (IF). The SISA should bound the true SISE with a certain confidence level. The IF is the quality indicator of the Signal-In-Space (SIS). If SISE cannot be bounded by SISA, the IF is set to “don’t use”. In a scalar representation of IF, the Signal-In-Space Monitoring Accuracy (SISMA) is also included. The SISMA quantity represents the accuracy of the estimation of the SISE.
Physically, the SISA and IF functions are implemented by two different chains, the navigation chain and integrity chain respectively (Dellago et al., Reference Dellago, Pieplu and Stalford2003). For the navigation chain, the observations collected at the Galileo Sensor Stations (GSS) are sent via a communication channel to the Orbit Synchronization Processing Facility (OSPF) at the Galileo Control Centre (GCC). The OSPF uses a sub-set of the data to compute the navigation message (including ephemeris and satellite clock data) and the SISA value for each of the satellites. Both the navigation message and the SISA are broadcast to users at a nominal rate of 100 minutes (i.e. the Age of Data (AoD) is kept at less than 100 minutes). This allows the contribution of the satellite ephemeris and clock errors to the User Equivalent Range Error (UERE) to be maintained at less than 65 cm (1σ) (Dellago et al., Reference Dellago, Pieplu and Stalford2003). For the integrity chain, another set of observations collected at the Galileo Sensor Station (GSS) are sent via a different communication channel to the Integrity Processing Facility (IPF) at the Galileo Control Centre (GCC). The IPF then uses the relevant data to estimate the SISE for each satellite and compares it with the corresponding SISA in real time. Based on the comparison, the IPF generates an integrity flag (including SISMA) and broadcasts it to users in real time. At the user side, both the received SISA and SISMA are used to calculate the Protection Level (as shown in Figure 1) which is compared against the alert limit.
Figure 1. Galileo integrity.
To protect all potential users, the SISE should be bounded by the SISA at the worst user location (WUL) otherwise an integrity alert should be raised. Note that the IF parameters broadcast to the user can contain any one of a ‘don’t use’ flag, a ‘not monitored’ flag or the SISMA. These are important to ensure both a high level of safety and service availability for all categories of users.
The SISE is a function of user location due mainly to the satellite orbit error. The location with the maximum error, WUL, must be found for each epoch to ensure that the IF covers all potential users. There are two schemes to cope with SISE; one is to find a SISE model on the surface of the Earth, which can fit the measured SIS errors at GSS stations, then to interpolate the SISE for any user location. The EGNOS (European Geostationary Navigation Overlay Service) type model and GIPA (Galileo Integrity Performance Assessment) models are based on this scheme (Werner et al., Reference Werner, Zink, Lohnert and Pielmeier2001). The surface could be a first order model (plane) or second order model (curved surface). The other SISE scheme is to calculate the satellite error from the measurements at the GSSs followed by the calculation of the SISE at the user location. Currently, WUL determination is based on a search on a defined grid. Depending on the SISE determination algorithm, the search space could be the footprint of the satellite with a certain grid density (e.g. EGNOS type model), or the border of satellite coverage at specific points for a given elevation angle (e.g. the GIPA first order model). The computation method (e.g. the GIPA second order model) is feasible only when the analytical maximum of SISE surface is within the visibility circle i.e. footprint (Werner et al., Reference Werner, Zink, Lohnert and Pielmeier2001).
The search method is a compromise between accuracy and search space (or computation time). A 0·5° search interval along the coverage border corresponding to circa 87 km is suggested by Werner et al., (Reference Werner, Zink, Lohnert and Pielmeier2001). A 1°×1° grid is used by Blomenhofer et al., (Reference Blomenhofer, Ehret, Leonard and Blomenhofer2004). Clearly, a smaller interval would be required to guarantee the accuracy of the determination of the WUL. This has the negative impact of increasing the number of points in the search space and hence the computation time. The Galileo system requires second by second computation of the IF. For global coverage with 27 satellites, computation time is clearly the limiting factor in the accuracy of the determination of the WUL given the strict requirements with respect to time-to-alert. Even if the GIPA second order model is used where the calculation of WUL is possible if it is within the satellite coverage, the maximum value of the SISE surface is not necessarily within the coverage border. In addition, the precise determination of the coverage border remains a problem for the border search method.
A new algorithm is proposed, which can calculate the WUL directly. The satellite related error estimation including the four-parameter and three-parameter models are analysed after the introduction. The proposed WUL algorithm is described in section 3. The results are presented and discussed in section 4. The paper is concluded in section 5.
2. SATELLITE RELATED ERROR ESTIMATION
Satellite errors can be expressed in either the four or three-parameter formats. The four-parameter format contains the three orbital components and the satellite clock error (Δx s Δy s Δz s Δt s) expressed in Earth-Centre-Earth-Fixed (ECEF) coordinates. The three-parameter format uses along-track (ΔL), cross-track (ΔC) and radial (ΔH) components, or along-track (ΔL), cross-track (ΔC) and clock (ΔT) components to express the satellite error.
In the four-parameter expression, the satellite error estimation approach is the same as the basic navigation method used for positioning, the difference being that the role of the satellites in the navigation solution is now played by the GSSs with the unknowns, the satellite orbit and clock errors (Δx s Δy s Δz s Δt s) expressed in the ECEF frame. A weighted least squares algorithm is used to determine the unknowns, and the corresponding quality and correlation parameters. The SISE and its quality parameter are then derived from these for each user-satellite pair.
The orbit error direction vector can be expressed as
where 
is the value of the satellite orbit error.
In the three-parameter expression, there are two ways to combine the radial and the clock error. One way is to merge the clock and the radial error as a radial error, and the other is to merge the clock and the radial error as a clock error. The SISE estimation procedure is the same as that for the four-parameter model but with one less unknown (either the clock error or the radical error). The SISMA represents the accuracy of the estimation of SISE at the WUL. It is derived from the accuracy estimates for each potential user.
3. WUL DETERMINATION
3.1. Philosophy of WUL determination
The satellite clock error is the same for all users within its coverage. Therefore the WUL depends only on the satellite orbit error. The idea behind the WUL computation is that the user at the worst location has the maximum range error. When the orbit error vector is extended, it intersects the Earth surface at the projection point. A user at the projection point has the maximum range error because the orbit error propagates to it completely as shown in Figure 2a. Hence, a user located at the projection point is in the worst location relative to the satellite orbit error. However, it is important to note that the projection point does not have to lie within the coverage of the satellite (footprint) as shown in Figure 2b. In this case (i.e. where the projection point is outside the footprint or even beyond the Earth) the WUL is on the boundary of the satellite footprint.
Figure 2a. Orbit error projection.
Figure 2b. Projection point outside footprint.
Under the assumption that the Earth is a sphere, the Earth centre, the satellite and the orbit error vector form a plane, which cuts the Earth as a great circle as shown in Figure 2c. The determination of the plane simplifies the problem from three dimensions (3D) to two dimensions (2D). The WUL is then at one point (or two points) on the great circle if the user is assumed to be on the surface of the Earth. Therefore, the WUL can be calculated by a reference point (known point) on the great circle, an azimuth (direction) from the reference point to WUL, and the distance between the reference point and WUL. The azimuth may be determined from two known points on the great circle. In this case, the Earth centre to the satellite vector can be used to determine one point on the great circle (foot point), and the satellite orbit error can be used to determine another point on the great circle. Either point can be considered to be a reference point. The two points can be used to determine the azimuth. The distance can be determined by the angle between the two vectors. The details are given in section 3.2. Clearly the WUL determined on the spherical model has to be transformed into the desired ellipsoidal model (section 3.4)
Figure 2c. The formation of great circle.
3.2. WUL computation for the four-parameter satellite error
3.2.1. WUL computation based on spherical model
The position of a satellite can be determined from the satellite’s almanac, broadcast ephemeris, or precise ephemeris. The Earth centre to the satellite direction vector can be expressed as
where x s, y s, z s are the position of satellite in ECEF coordinates, is the distance between the centre of the Earth and the satellite.
For the spherical model, the location of the satellite foot point can be calculated from
where, φ f, λ f are the latitude and longitude of the satellite foot point respectively.
The basic relationships between satellite coverage and various angles are shown in Figure 3. The nadir, central and elevation angles are denoted by α, β, and θ respectively. r e is the radius of the Earth. α 0 is the nadir angle for a given mask angle θ 0, β 0 is the corresponding central angle for the same mask angle. The central angle is the distance in radians between the foot point F and the worst user location W.
Figure 3. Relationships between satellite coverage and various angles.
The nadir angle α can be calculated based on the satellite orbit error described in section 2. In order to do so, the satellite orbit error direction vector is shifted to the centre of the Earth. Therefore, the location of the point A(or A′) expressed as φ A, λ A in terms of latitude and longitude can be found on the surface of the Earth corresponding to error vector SE (or SE′) by equations (3) and (4) from the value of orbit error direction vector e E. γ is the arc cosine of the dot product of the centre of the Earth to the satellite direction vector e S and the satellite orbit error direction vector e E (γ takes the values between 0 and 180° inclusive). α has the same value as γ when γ is less than 90°, or the value of supplementary of γ when γ is larger than 90°.
The fundamental relationships between these angles and distances are
The radius of the Earth r e is known. Furthermore, the distance between the satellite and the Earth’s centre r sat can be calculated.
The azimuth of can be calculated from the location of points A and F by
Under the condition that α⩽α 0, the sum of α, β and θ is 90°. The value of any angle can be determined from the values of the other two angles. Under the condition that α>α 0 the projection point is outside the satellite coverage, and β>β 0. In this case, the worst user location is on the border of the satellite coverage i.e. β=β 0.
The WUL point W can be calculated by using:
where φ W, λ W are the latitude and longitude of point W and
3.2.2. Azimuth correction
The direction vector of the WUL on the spherical model can be expressed as
Where, , (X wY wZ w) is the WUL in ECEF coordinates.
Theoretically, the direction vectors e WS, e S and e E should be in one plane. However, due to the error in the computation of direction , φ W and λ W, these direction vectors may not be in one plane, which results in WUL errors even in the spherical model and causes a further problem when transformed to the ellipsoidal model.
The error in the azimuth (expression 8) must be corrected to ensure that these three vectors are in one plane before transforming the WUL to the ellipsoidal model.
The direction vector perpendicular to the e WS, e S plane can be determined from:
Where×denotes cross product, |e WS×e S| denotes the norm of e WS×e S.
Similarly, the direction vector e SE perpendicular to the e S, e E plane can be determined.
Therefore, the error in direction can be determined by:
Where, ● denotes dot product,
By adding the correction to and feeding it back to expression (9) and (10), a new WUL in the spherical model can be determined and should be in the plane determined by e S and e E.
3.3. WUL computation for the three-parameter satellite error
In some cases, such as the absence of sufficient measurements, three-parameter satellite error models are introduced by combining the orbit radial and the clock errors. There are two WUL computation algorithms corresponding to the two ways of combining the radial and clock errors (section 2).
In the case where the satellite error vector is expressed in terms of satellite along-track (ΔL), cross-track (ΔC) and radial (ΔH) components, the same scheme described in section 3.2 is applied to compute the WUL when the three-parameter expression is transformed to ECEF coordinates.
In the case where the satellite error vector is expressed in terms of satellite along-track (ΔL), cross-track (ΔC) and clock (ΔT) components, as mentioned above, the satellite clock error is the same for all users within its coverage. Hence, the orbit error here refers only to the orbit cross track and along track error components. The orbit error vector is almost perpendicular to the centre of the Earth to satellite vector (i.e. as shown γ≈90°), because the satellite orbit is roughly a circle due to the small value of the orbit eccentricity. In this case, the sign of cos(γ) in expressions (9) and (10) takes both positive and negative forms; one represents the maximum positive error value, and the other the maximum negative error value. There are two WULs on the Earth surface, both lie on the boundary of the satellite footprint. The WULs determination scheme is the same as described in section 3.2 with the central angle to be β 0.
3.4. Transformation to ellipsoidal model
The WUL determined in the spherical model is denoted as WULS in Figure 4. Any point that lies on the line determined by the WULS and the satellite is a potential WUL. The difference is the distance between the point and the centre of the Earth, which also reflects the height of user. The user is assumed to be on the surface of the Earth. The WUL point can be expressed in either the spherical or the ellipsoidal model. However, the ellipsoidal model expression is required here because it is a more accurate approximation of the Earth. The transformation is performed by an estimation of the WUL followed by a distance correction if necessary.
Figure 4. Transformation to ellipsoidal model.
3.4.1. WUL estimation in ellipsoidal model
The radius of curvature in the prime vertical can be estimated as:
where R e is the semi-major axis, f is the flattening, φ w is the latitude in the spherical model and used to estimate the initial value of R N.
The WUL Estimation in Figure 4 is determined by extending the line joining the centre of the Earth to the WULS to the point where it intersects the surface of the ellipsoidal model of the Earth. Therefore, the distance between the centre of the Earth and WUL Estimation can be calculated from
![R_{OW} \equals R_{N} \left\vert {\left[ {cos\phi _{w} cos\lambda _{w} \quad cos\phi _{w} sin\lambda _{w} \quad \lpar 1 \minus f\rpar ^{\setnum{2}} sin\phi _{w} } \right]} \right\vert.](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160913095525830-0662:S037346330600381X_eqn15.gif?pub-status=live){kind=small}
R OW is fed back into expression (5) to replace r e and the central angle is adjusted to fit the value of R OW. Therefore, the WUL Estimation on the spherical model with radius as R OW can be calculated, which is an estimation of WUL in the ellipsoidal model since it is quite close to the theoretical WUL.
The location of the WUL can then be transformed to ECEF coordinates by
![\left[ {x_{w} \quad y_{w} \quad z_{w} } \right]^{T} \equals \left[ {R_{OW} cos\phi _{w} cos\lambda _{w} \quad R_{OW} cos\phi _{w} sin\lambda _{w} \quad R_{OW} sin\phi _{w} } \right]^{T}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160913095525830-0662:S037346330600381X_eqn16.gif?pub-status=live){kind=small}
Where (x w, y w, z w) is the position of WUL in ECEF coordinates which can be further converted to the geodetic location in the ellipsoidal model by normal coordinate conversion methods.
3.4.2. Distance (Central Angle) Correction
The WULEstimation in ECEF coordinates estimated in section 3.4.1 may contain a small error. It can be checked by comparing the WUL to the satellite vector e WULS and the satellite orbit error vector e E. The angular error Δα is the angle between these two vectors as shown in Figure 4. It is the arc cosine of the dot product of the vector e WULS and e E. The equivalent error in distance is roughly the product of angular error Δα and the distance between the satellite and estimated WUL. If the distance error is larger than the required accuracy e.g. 100 metres, Δα can then be used to correct the central angle β and then fed back to expressions (9) and (10). Iteration may be needed until the required accuracy is met, leading to the final WUL. The distance correction is carried out in the plane determined in section 3.2.2. The process does not introduce any direction errors; there is no need for a direction correction in the distance correction.
3.5. Boundary determination
If the projection point is outside the satellite footprint initially, the WULs is on the boundary of the footprint. In this case, the error Δα can no longer be used to correct the distance. Instead, the error between the elevation angle at WULEstimation in the ellipsoidal model and mask angle can be used to adjust the central angle β and fed back to expressions (9) and (10).
Problems may arise when the WUL determined using the method described in section 3.3 is transformed to the ellipsoidal model if the projection point is near the footprint boundary. Considering that α 0 in Figure 3 is given by the mask angle θ 0 in the spherical model, the transformation may result in a shift of the projection point from the inside of the spherical footprint to the outside of the ellipsoidal footprint, or vice versa. The final result is that an amount of error remains in the calculated WUL. Therefore, the elevation angle must also be compared to the mask angle before the correction. If the elevation angle is larger than the mask angle, the case that the projection point shifts from the outside of the footprint to the inside occurs. Therefore, the projection point is treated as inside the satellite footprint.
On the other hand, if the projection point is within the satellite footprint initially, the WUL can be calculated using the above method. However, more steps are needed. The computation of the elevation angle at the location of the projection point in the ellipsoidal model is required. If the elevation angle is larger than the mask angle, the location of the projection point is the WUL. If the elevation angle is smaller than the mask angle, the case that the projection point shifts from the inside of the footprint to the outside occurs. Therefore, the WUL is on the boundary of the footprint and can be calculated using the method described in the paragraph above.
3.6. Algorithm Implementation Steps
The implementation steps of the algorithm are presented in Figure 5.
Figure 5. Algorithm implementation of WUL computation.
4. RESULTS
To have a full understanding of the WUL, the range error distribution corresponding to the unit error in satellite orbit has been analysed. Furthermore, a number of case studies have been carried out to verify the performance of the new WUL algorithm. Both considered all possible cases including:
● the projection point is exactly at the satellite foot point,
● the projection point is within the satellite footprint
● the projection points are exactly on the border of the satellite footprint
● the projection point is outside the satellite footprint including points beyond the Earth.
4.1. Range Error distribution
Four scenarios were investigated to characterise the range error distributions corresponding to the unit error in satellite orbit. The angle γ defined in section 3 is the key parameter for analysing these four scenarios.
● γ=0, where the orbit error unit vector is the same as the unit vector of the centre of the Earth to the satellite. The error distribution is shown in Figure 6a. The user location with the maximum error is at the satellite foot point, and degrades with increasing distance from the satellite foot point.
● γ<α 0. The user location with the maximum error is within the coverage of the satellite (footprint). The unit error distribution is shown in Figure 6b. The peak in the Figure is the WUL.
● γ=90°. This is the same as the case when satellite orbit radial error is absorbed into the clock error. The error distribution is shown in Figure 6c. Two peaks occur on the satellite coverage border, one with a maximum positive value, and the other with a maximum negative value.
● 90°<γ<(180−α 0). The error distribution is shown in Figure 6d. The shape of the figure is similar to Figure 6c. However, the error value is different. In this case, only one peak occurs on the satellite coverage border.
Figure 6a: Error distribution case 1.
Figure 6b: Error distribution case 2.
Figure 6c: Error distribution case 3.
Figure: 6d Error distribution case 4.
4.2. Case studies
The following four cases were used to investigate the performance of the proposed algorithm:
● The orbit error vector is the same as the vector of the centre of the Earth to the satellite.
● The projection of the orbit error vector is within the satellite coverage of 15° mask angle.
● The orbit error vector is perpendicular to the vector of the centre of the Earth to the satellite.
● The projection of the orbit error vector is outside the satellite coverage of 15° mask angle.
Table 1 shows the results of the computed WUL, distance error, elevation angle, and the number of iterations for the assumed the error vector and a satellite position 1 at (7981408.9299816 16578447.5771181 23186465.6087373) in ECEF coordinates.
Table 1. Results for satellite position 1.
Table 2 shows the results of the computed WUL, distance error, elevation angle, and the relevant number of iterations for the assumed the error vector and a satellite position 2 at (29556256.5793380 1607884.0619139 49049.5645304) in ECEF coordinates.
Table 2. Results for satellite position 2.
For the cases where the WUL is within the satellite footprint, the accuracy of the proposed method can reach the metre level with only one iteration in the correction of the direction and one iteration in the correction of the distance. However, in the current grid search-based method, since the central angle is roughly 63° for a 15°mask angle, the search method requires as much as 15876 sample points for a 1°×1° grid, while the great circle distance error could be as large as 78 km. For the cases where the WUL is outside the satellite footprint, the accuracy of the proposed method can reach the metre level with one iteration in the correction of the direction and a few iterations in the correction of the distance. On the other hand, the search-based method requires 360 sample points for one degree interval around the boundary, while the great circle distance error could reach 55 km. Of course the accuracy of the search-based method can be improved by reducing the sample (or grid) interval. However, as has been explained in the paper, their accuracy is limited by the computational workload required to deal with a more granular grid.
The case studies demonstrate the accuracy and efficiency of the algorithm.
5. CONCLUSIONS
The proposed algorithm uses the resolved orbital errors to calculate an approximate WUL on a spherical model of the Earth. The next step is to determine the WUL on a specified ellipsoid model using an iterative transformation process. The method is more accurate and computationally efficient than the current grid search-based methods. The accuracy is achieved through the correction of the direction and distance. The computational workload is reduced by avoiding extensive searching. Both three-parameter and four-parameter expressions of the satellite errors are considered. The possible overpass boundary problem caused by the transformation from the spherical model to the ellipsoidal model is also addressed. The advantages (accuracy and temporal efficiency) make it valuable in meeting the real-time requirements for the generation of the Galileo Integrity Flag.
