Covariance Analysis of Real-Time Precise GPS Orbit Estimated from Double-Di ﬀ erenced Carrier Phase Observations

: The covariance of real-time global positioning system (GPS) orbits has been drawing attention in various ﬁelds such as user integrity, navigation performance improvement, and fault detection. The international global navigation satellite system (GNSS) service (IGS) provides real-time orbit standard deviations without correlations between the axes. However, without correlation information, the provided covariance cannot assure the performance of the orbit product, which would, in turn, causes signiﬁcant problems in fault detection and user integrity. Therefore, we studied real-time GPS orbit covariance characteristics along various coordinates to e ﬀ ectively provide conservative covariance. To this end, the covariance and precise orbits are estimated by means of an extended Kalman ﬁlter using double-di ﬀ erenced carrier phase observations of 61 IGS reference stations. Furthermore, we propose a new method for providing covariance to minimize loss of correlation. The method adopted by the IGS, which neglects correlation, requires 4.5 times the size of the covariance to bind orbit errors. By comparison, our proposed method reduces this size from 4.5 to 1.3 using only one additional parameter. In conclusion, the proposed method e ﬀ ectively provides covariance to users. all covariance ellipsoids of the RSW and RAC appear more uniform than those of ECI and ECEF. The estimated covariance in the RSW and RAC coordinate systems appear in the form of a long ellipsoid in the along-track direction. For the remainder of the study, we focused on the analysis to RSW and RAC frames.


Introduction
In recent years, precise real-time global positioning system (GPS) orbit products have extended the capability to support real-time applications such as autonomous driving, intelligent transportation systems, and collision avoidance [1,2]. The real-time precision of orbits enables navigation systems to overcome problems arising from orbit errors in real-time observations [3].
Precise real-time GPS orbits have been determined with the international global navigation satellite system (GNSS) service (IGS) forming the mainstream service. IGS provides ultra-rapid orbits as real-time precise orbits [4]. Ultra-rapid orbits are determined using recent satellite arcs of 3 days to predict the orbit for 24 h thereafter [5]. The maximum accuracy of the predicted orbit is 5 cm at 1D root mean square (RMS). Another such service is the IGS real-time service, which is supported by ten analysis centers [2]. One of the analysis centers, the Centre National d'Études Spatiales (CNES), estimates the satellite orbit and clock together using undifferenced GPS observations. Other IGS analysis centers often use ultra-rapid solutions for precise orbits and concentrate on precisely estimating satellite clocks [6][7][8]. Furthermore, precise point positioning (PPP)-based commercial services have similar strategies to generate precise orbits and clocks when compared with IGS real-time products. Past studies have focused on verifying the accuracy of these real-time precise orbits [5,8,9].
Interest in the covariance of real-time precise GPS orbits, as well as in the accuracy, has increased in the context of user integrity, navigation performance improvement, and fault detection [10][11][12][13]. In terms of safety, covariance is one of the most important factors that provide integrity that ensures correct position information. Range error due to orbit error should be overbounded to ensure user position integrity [14,15]. In general, systems provide the accuracy of signal-in-space range error (SISRE) or orbit full covariance. The accuracy of SISRE is calculated with proper weightings for radial, along-track, and cross-track standard deviations [16,17]. Providing full covariance enables each user to propagate the accuracy of SISRE by projecting the error ellipsoid along its line of sight [11], which can reduce specific accuracy of SISRE according to the various positions of the users [12]. Therefore, to generate appropriate accuracy of SISRE, proper error ellipsoid which contains the true orbit should be obtained.
However, orbit full covariance is not provided from real-time orbit products. Previous studies utilized the maximum accuracy or the stochastics of orbit errors to analyze the effect of covariance usage. El-Mowafy suggested using the orbit covariance to detect faults or meaconing errors in IGS RTS correction and demonstrated the significant advantage of a new fault detection model over traditional models under a meaconing attack [18]; the expected effect of covariance was considered using the maximum accuracy of each axis, which is greater than the full covariance. The full covariance would provide orbit uncertainty more appropriately and would improve the performance of fault detection. In addition, Cheng et al. [19] analyzed the user range accuracy performance of real-time ephemeris. They studied the characteristics of long-term error stochastics over a year to provide the performance of user range accuracy. However, as only the average characteristics can be obtained over a long duration, the correlation for real-time covariance should be analyzed for real-time applications. Therefore, we estimated the real-time orbit covariance to identify more realistic covariance characteristics.
To support covariance-based applications, IGS plans to provide real-time full covariance, although it is not expected to be provided shortly in the near future [20]. The real-time orbit standard deviation is presently provided over the XYZ components of the Earth-centered Earth-fixed (ECEF) frame, without any correlations between the axes. However, if the correlation information is neglected, the provided accuracy cannot ensure the performance of the orbit product, which in turn causes significant problems in fault detection and user integrity. Therefore, it is safe to provide a conservative representation of the orbit error distribution, however unnecessary overbounding decreases the availability of the system. This paper proposes an effective covariance provision method considering the correlations of real-time GPS orbits. We analyzed these correlations using a real-time GPS precise orbit estimator and studied the covariance along various coordinates. The results demonstrate the real-time characteristics of the correlations, which cannot be determined based on long-term analyses. Considering the correlation, we propose a covariance provision method and evaluate it using the number of parameters and ratio of the provided covariance volume to the full covariance volume.
The remainder of this paper is organized as follows. Section 2 presents the details of the orbit determination tool and the theoretic background of real time orbit covariance due to orbit dynamics; Section 3 discusses the experimental results; Section 3.1 verifies the orbit determination system in relation to the IGS final orbit; Section 3.2 presents the experimental results of covariance analysis and proposes a new frame to minimize correlations between the axes; Section 3.3 assesses several covariance provision methods; Section 4 discusses our findings; and finally, Section 5 presents our conclusions.

Observational Data
This study uses double-differenced carrier phase (DDCP) to focus on orbit determination without clock errors. To obtain cm-level accuracy, CP observations are generally used to determine precise orbits. CP is defined as the phase difference between the satellite and the receiver, and it contains Remote Sens. 2019, 11, 2271 3 of 18 satellite clock errors and receiver clock errors in addition to the distance information between the satellite and receiver [21]. Additionally, there are phase center offsets (PCO) and variations (PCV), phase wind up (PWU) and tidal effects, which are generally corrected by models for GPS precise orbit determination [22]. Researchers normally use these observations in the form of un-differenced CP (UDCP) or DDCP data. Clock errors can be eliminated by means of the double-differencing method. Furthermore, if the network comprises the same type of receiver, the hardware bias is also eliminated with the double-differencing method. In contrast, the noise of the receiver is doubled. Table 1 compares the error components of UDCP and DDCP, wherein 'V' indicates the presence of an error component.
set to 5 • . The IERS-predicted EOP is used for coordination transformation, and broadcast ephemeris is used for preprocessing.
filter are raw measurements, the broadcast (BRDC) orbit, and the International Earth Rotation and Reference System (IERS)-predicted earth orientation parameters (EOP). Raw measurements are acquired in the RINEX format from the IGS data archive [22]. Dual-frequency PR and CP data are acquired from 61 globally distributed stations across 30 s intervals. Figure 2 shows the station locations and baselines. The stations are chosen to form a station network with homogeneous receivers, Trimble receiver in this case, in order to avoid estimation of hardware bias errors. The mask angle is set to 5°. The IERS-predicted EOP is used for coordination transformation, and broadcast ephemeris is used for preprocessing.  receivers, Trimble receiver in this case, in order to avoid estimation of hardware bias errors. The mask angle is set to 5°. The IERS-predicted EOP is used for coordination transformation, and broadcast ephemeris is used for preprocessing.  First, the raw data are preprocessed to detect and compensate for cycle slips. Subsequently, the main filter uses the 'cleaned' CP to determine the orbit and its covariance. The main filter state consists of the satellite position, vector, coefficients of solar radiation pressure, wet zenith tropospheric delays, and ambiguities. For more accurate orbit results, the main filter should resolve the L2 ambiguities. In this step, the fixed WL ambiguity is required. Therefore, the WL ambiguity filter estimates the WL ambiguity using DDCPMW observations beforehand. To mitigate the effect of multipath, this filter also estimates the multipath error as well as WL ambiguity by modeling the first-order Gauss-Markov process [24,25]. The detailed method is described in [24]. Finally, the GPS orbits and their covariance are generated in real-time. Tables 2 and 3 show the states of the WL ambiguity filter and main Kalman filter.  In order to generate precise orbits, the cycle slip should be detected and compensated. Cycle slip is a discontinuity that appears when a receiver experiences a loss of lock [21]. This study adopted the preprocessing algorithm proposed in the Bernese GNSS software [26]. First, a DDIF observation was fitted with the second-order polynomial for outlier and cycle slips. Second, the time-differenced MW and ionosphere combinations were utilized to detect and compensate for the cycle slip. The algorithm was verified in reference [24].
Several additional corrections are considered in this preprocessing step. When the baseline length is small, errors such as phase wind up and tidal displacements can be eliminated with the double differencing method [3,27]. However, the baselines of our network, which range from 120 to 3800 km, are too long to neglect the extra errors. Therefore, we compensated the errors; Table 4 summarizes the corresponding strategy [28][29][30][31]. The strategies use the GPS attitude model for nominal attitude [28], which is also called yaw-steering attitude.

Main Filter Description
The main filter utilizes a sequential extended Kalman filter (EKF) for real-time orbit determination. The sequential approach processes incoming data immediately, and can quickly respond to and reflect the real-time orbit [32]. The EKF utilizes two steps: time update (TU) and measurement update (MU) steps. The orbit and covariance are estimated over repeated TU and MU iterations. The orbit is initialized with the broadcast orbit, and the initial standard deviation is set to 2 m for each axis. In the TU step, the previous step state is propagated with system dynamics. The state of the satellite position and velocity are propagated by orbit dynamics. On the other hand, the troposphere wet zenith delay of every station makes use of the first-order Gauss-Markov process as the propagation model. The coefficients of solar radiation pressure and L2 integer ambiguities are estimated together. The state transition matrix of the orbit position and velocity requires numerical derivatives. The setting values of filter is described in Tables 5 and 6.  The MU step is performed using the DDCPIF observations of the network. To estimate the troposphere wet zenith delay, we eliminate the dry component using the Saastamoinen model [33]. Furthermore, the dry and wet mapping function is based on the Niell mapping function [34].

Orbit Propagation
Each satellite orbit is propagated by means of the Runge-Kutta 68 numerical integral method [32] using satellite dynamics. The integration step is 30 s. For the orbit propagation, the state and dynamics utilize the J2000 coordinate system. The equinox-based transformation is implemented using IAU-2000A as per [35]. The satellite acceleration is calculated with the use of Cowell's method [36], which adds perturbing accelerations to the two-body equation to generate a more accurate equation of motion. The satellite acceleration is expressed, as follows: where → r and → a the position and acceleration, respectively. We add the perturbing effects of the non-spherical central body, the third-body effects of the Sun and Moon, solar radiation pressure, and the tides. The Earth's gravitation field is implemented as per the earth gravity model (EGM) 96 and is truncated to the 12th degree and order for GPS orbits [37]. Third-body effects consider the Sun and Moon as the third planet, and utilizes the DE405 planet model [38]. Tidal and general relativistic effects are examined following IERS 2010 standards. For solar radiation pressure, precise GPS orbit products estimate the coefficients of the model of dynamics in general. The Empirical Center for Orbit Determination in the Europe Orbit Model (ECOM) [39] is a popular dynamics model for real-time GPS orbit products [40,41]. We applied the ECOM model with nine radiation pressure terms for each satellite. The nine coefficients per satellite are estimated together. In addition, the shadows of the Earth and Moon comprise a conical model [42]. Table 7 summarizes the dynamic model.

Propagation Characteristics of Orbit Dynamics
Before the experimental analysis, we briefly explain the characteristics of covariance theoretically by means of orbit dynamics. In general, the orbit error can be induced as the relative movements of two objects, where one orbit is on the true orbit and the other is on the estimated orbit. The orbit propagation properties have been studied using relative motion equations in the RSW frame [43][44][45].
The geometric modeling of relative motion can express orbit error by means of orbital elements as follows for a near-circular orbit [45]: x ≈ δa − acos(f)δe (7) y where, x, y, z, v x , v y , and v z represent the orbit errors corresponding to the position and velocity vector for each axis in the RSW frame; a, e, i, Ω, w, f and M denote the orbital elements of the semi-major axis, eccentricity, inclination angle, right ascension of the ascending node, argument of perigee, argument of latitude, true anomaly, and mean anomaly, respectively; n denotes the mean motion; and δ is the error of each component. The semi-major axis error δa is affected by disturbances, and it is generally modeled by a linear function of time.
It affects mean anomaly M, which is expressed as: Upon utilizing a linear function of, δa, x, y, z, v x , and v y can be, respectively, expressed as: where τ denotes the elapsed time with propagation, and subscript 0 indicates the initial value; furthermore, δ . a denotes the rate of the semi-major axis error. Equations can aid us in understanding the orbit propagation characteristics. Errors along the radial and along-track directions consist of secular, periodic, and constant components. The square of parameter τ diverges rapidly as time elapses, and this leads to large variance of the along-track errors. The periodic terms are due to δe and their period is the same as the orbit period. In general, δe, δi, δΩ, and δw are all extremely small.
In addition, the correlation between the radial and along-track values can be explained using Equation (19). The secular component of the along-track error exhibits a proportional relationship with the radial velocity error, and the proportional constant has a negative value and is expressed as: This leads to negative correlation of the errors of the radial and along-track directions. The cross-track position and velocity errors are periodic for the argument of latitude u.

Verification of Orbit Filter and Covariance Analysis
First, we present our evaluation of the position results in relation to the IGS final orbits (Section 3.1). Section 3.2 presents the error covariance characteristics over time and for several frames; IGS final orbits are not used. Furthermore, the error correlations between the axes are analyzed. In Section 3.3, these characteristics are utilized to propose and evaluate a new covariance parameterizing method.

Experimental Environment
We obtained the results of the orbit filter using the dual-frequency PR and CP of the 61 IGS stations shown in Figure 2; these were collected on 15-16 January 2018. The filter generated real-time orbits at intervals of 30 s. In this section, we present our verification of the filter performance in relation to the IGS final product. The accuracy of IGS final orbit products [46], which are the most precise GPS orbits, is 2.5 cm [47]. The position error is the difference between the estimated satellite positions and the IGS final product. In addition, a simulation was conducted to verify the filter. The simulation data are summarized in Table 8.  Figure 3 shows the time history of the pseudo-random noise (PRN) 9 orbit error of simulation and real data. The red, green, and blue curves denote the 3D position error, radial error, and estimated standard deviation of the 3D position error, respectively. The initial orbit conditions are calculated using the BRDC orbits, and therefore they exhibit meter-level accuracy. The orbit accuracy convergences to the centimeter level after 24 h. Figure 3 shows the time history of the pseudo-random noise (PRN) 9 orbit error of simulation and real data. The red, green, and blue curves denote the 3D position error, radial error, and estimated standard deviation of the 3D position error, respectively. The initial orbit conditions are calculated using the BRDC orbits, and therefore they exhibit meter-level accuracy. The orbit accuracy convergences to the centimeter level after 24 h.

Actual Error Distribution of Real-Time Precise Orbit
We analyzed the actual error distribution of real-time precise orbits to confirm that the error distribution is properly bounded to the estimated covariance. Normalized error, defined as the value of the residual divided by the estimated standard deviation, are analyzed over a 48 h interval to indicate the error distribution. We performed not only the probability density function (PDF) bounding analysis but also the cumulative distribution function (CDF) bounding analysis, a method used to show how well the tail distribution is bounded [40,41].
The error distribution of IGS ultra-rapid was also analyzed and compared. IGS ultra-rapid provides real-time covariance information without correlations. Figure 5 shows that the error distribution is not bounded by the covariance information. Thus, the application of covariance will suffer from unbounded orbit errors. Owing to this, users cannot confidently rely on the covariance. To prevent covariance-based applications from failing, users should apply a conservative error distribution using a scaling factor.  Figure 4 shows the root mean square (RMS) orbit errors of all satellites during the last 24 h of a given day. The mean 3D RMS value is 7.8 cm, and the radial track error is~2 cm. The along-track error was larger than radial and cross-track errors. Figure 3 shows the time history of the pseudo-random noise (PRN) 9 orbit error of simulation and real data. The red, green, and blue curves denote the 3D position error, radial error, and estimated standard deviation of the 3D position error, respectively. The initial orbit conditions are calculated using the BRDC orbits, and therefore they exhibit meter-level accuracy. The orbit accuracy convergences to the centimeter level after 24 h.

Actual Error Distribution of Real-Time Precise Orbit
We analyzed the actual error distribution of real-time precise orbits to confirm that the error distribution is properly bounded to the estimated covariance. Normalized error, defined as the value of the residual divided by the estimated standard deviation, are analyzed over a 48 h interval to indicate the error distribution. We performed not only the probability density function (PDF) bounding analysis but also the cumulative distribution function (CDF) bounding analysis, a method used to show how well the tail distribution is bounded [40,41].
The error distribution of IGS ultra-rapid was also analyzed and compared. IGS ultra-rapid provides real-time covariance information without correlations. Figure 5 shows that the error distribution is not bounded by the covariance information. Thus, the application of covariance will suffer from unbounded orbit errors. Owing to this, users cannot confidently rely on the covariance. To prevent covariance-based applications from failing, users should apply a conservative error distribution using a scaling factor.

Actual Error Distribution of Real-Time Precise Orbit
We analyzed the actual error distribution of real-time precise orbits to confirm that the error distribution is properly bounded to the estimated covariance. Normalized error, defined as the value of the residual divided by the estimated standard deviation, are analyzed over a 48 h interval to indicate the error distribution. We performed not only the probability density function (PDF) bounding analysis but also the cumulative distribution function (CDF) bounding analysis, a method used to show how well the tail distribution is bounded [40,41].
The error distribution of IGS ultra-rapid was also analyzed and compared. IGS ultra-rapid provides real-time covariance information without correlations. Figure 5 shows that the error distribution is not bounded by the covariance information. Thus, the application of covariance will suffer from unbounded orbit errors. Owing to this, users cannot confidently rely on the covariance. To prevent covariance-based applications from failing, users should apply a conservative error distribution using a scaling factor.
In Figures 6 and 7, panels (a) to (c) correspond to the X-, Y-, and Z-axis results of the Earth-centered inertial (ECI) frame. Figure 6 shows the PDF of the normalized error. The blue and green bars indicate the actual sample simulation PDF data and real data, respectively, while the red curve indicates a normal distribution. Figure 7 represents the CDF for each axis. The actual simulation CDF data and real data are indicated by the blue and green lines, respectively. The red curve represents the normal distribution curve. The normal distribution of the covariance guarantees the conservative distribution In Figure 6 and Figure 7, panels (a) to (c) correspond to the X-, Y-, and Z-axis results of the Earth-centered inertial (ECI) frame. Figure 6 shows the PDF of the normalized error. The blue and green bars indicate the actual sample simulation PDF data and real data, respectively, while the red curve indicates a normal distribution. Figure 7 represents the CDF for each axis. The actual simulation CDF data and real data are indicated by the blue and green lines, respectively. The red curve represents the normal distribution curve. The normal distribution of the covariance guarantees the conservative distribution of orbit errors for each axis. Figure 6 and Figure 7 confirm that the error distribution is properly bounded to the estimated covariance.

Orbit Covariance in Several Coordinates
The covariance characteristics of satellite orbits were analyzed to determine the best coordinate system to minimize the correlation between each axis. We studied the covariance in four different frames: ECI, ECEF, radial-transverse-normal (RSW), and radial along-and cross-track (RAC). The RSW frame is a satellite-fixed coordinate system that is defined using the satellite position and velocity vector. The R-axis represents the radial track, which refers to the direction towards the satellite from the center of the Earth. The W-axis represents the cross-track, which is perpendicular to both the satellite position vector and the velocity vector. The S-axis is defined as the along-track,  In Figure 6 and Figure 7, panels (a) to (c) correspond to the X-, Y-, and Z-axis results of the Earth-centered inertial (ECI) frame. Figure 6 shows the PDF of the normalized error. The blue and green bars indicate the actual sample simulation PDF data and real data, respectively, while the red curve indicates a normal distribution. Figure 7 represents the CDF for each axis. The actual simulation CDF data and real data are indicated by the blue and green lines, respectively. The red curve represents the normal distribution curve. The normal distribution of the covariance guarantees the conservative distribution of orbit errors for each axis. Figure 6 and Figure 7 confirm that the error distribution is properly bounded to the estimated covariance.

Orbit Covariance in Several Coordinates
The covariance characteristics of satellite orbits were analyzed to determine the best coordinate system to minimize the correlation between each axis. We studied the covariance in four different frames: ECI, ECEF, radial-transverse-normal (RSW), and radial along-and cross-track (RAC). The RSW frame is a satellite-fixed coordinate system that is defined using the satellite position and velocity vector. The R-axis represents the radial track, which refers to the direction towards the satellite from the center of the Earth. The W-axis represents the cross-track, which is perpendicular to both the satellite position vector and the velocity vector. The S-axis is defined as the along-track,  In Figure 6 and Figure 7, panels (a) to (c) correspond to the X-, Y-, and Z-axis results of the Earth-centered inertial (ECI) frame. Figure 6 shows the PDF of the normalized error. The blue and green bars indicate the actual sample simulation PDF data and real data, respectively, while the red curve indicates a normal distribution. Figure 7 represents the CDF for each axis. The actual simulation CDF data and real data are indicated by the blue and green lines, respectively. The red curve represents the normal distribution curve. The normal distribution of the covariance guarantees the conservative distribution of orbit errors for each axis. Figure 6 and Figure 7 confirm that the error distribution is properly bounded to the estimated covariance.

Orbit Covariance in Several Coordinates
The covariance characteristics of satellite orbits were analyzed to determine the best coordinate system to minimize the correlation between each axis. We studied the covariance in four different frames: ECI, ECEF, radial-transverse-normal (RSW), and radial along-and cross-track (RAC). The RSW frame is a satellite-fixed coordinate system that is defined using the satellite position and velocity vector. The R-axis represents the radial track, which refers to the direction towards the satellite from the center of the Earth. The W-axis represents the cross-track, which is perpendicular to both the satellite position vector and the velocity vector. The S-axis is defined as the along-track,

Orbit Covariance in Several Coordinates
The covariance characteristics of satellite orbits were analyzed to determine the best coordinate system to minimize the correlation between each axis. We studied the covariance in four different frames: ECI, ECEF, radial-transverse-normal (RSW), and radial along-and cross-track (RAC). The RSW frame is a satellite-fixed coordinate system that is defined using the satellite position and velocity vector. The R-axis represents the radial track, which refers to the direction towards the satellite from the center of the Earth. The W-axis represents the cross-track, which is perpendicular to both the satellite position vector and the velocity vector. The S-axis is defined as the along-track, and it is perpendicular to the R-and W-axes. Given the satellite position and velocity vector, the RSW frame can be defined as follows [49]: where → r ECI and → v ECI represent the position and velocity in the ECI frame, respectively. In general, the RSW frame is calculated using the ECI frame states. However, a user only has a satellite's position and velocity vector as per the ECEF frame. To generate the RSW frame using ECI vectors, ECEF should be transformed to ECI. To this end, users are required to overcome the associated computational burden; this is a substantial problem, particularly for low-cost user systems such as smart phones. Furthermore, additional communication is required to obtain EOP parameters. Therefore, IGS RTS products are provided in an RAC frame, and these are calculated utilizing the ECEF position and velocity vectors. The RAC frame is similar to the RSW frame based on ECEF states. We define the RAC frame in the following section. In this context, we note that in this study, we examined the covariance characteristics using not only the general RSW frame, but also an RAC frame calculated using the ECEF states: where → r ECEF and → v ECEF denote the position and velocity, respectively, in the ECEF frame. Figure 8 shows the entire satellite covariance after 23 h in the four different frames (ECI, ECEF, RSW, and RAC). In Figure 8a,b, the axis-to-axis correlation of a satellite is different from that of other satellites. On the other hand, all covariance ellipsoids of the RSW and RAC appear more uniform than those of ECI and ECEF. The estimated covariance in the RSW and RAC coordinate systems appear in the form of a long ellipsoid in the along-track direction. For the remainder of the study, we focused on the analysis to RSW and RAC frames.
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 18 and it is perpendicular to the R-and W-axes. Given the satellite position and velocity vector, the RSW frame can be defined as follows [49]: where ⃗ and ⃗ represent the position and velocity in the ECI frame, respectively. In general, the RSW frame is calculated using the ECI frame states. However, a user only has a satellite's position and velocity vector as per the ECEF frame. To generate the RSW frame using ECI vectors, ECEF should be transformed to ECI. To this end, users are required to overcome the associated computational burden; this is a substantial problem, particularly for low-cost user systems such as smart phones. Furthermore, additional communication is required to obtain EOP parameters. Therefore, IGS RTS products are provided in an RAC frame, and these are calculated utilizing the ECEF position and velocity vectors. The RAC frame is similar to the RSW frame based on ECEF states. We define the RAC frame in the following section. In this context, we note that in this study, we examined the covariance characteristics using not only the general RSW frame, but also an RAC frame calculated using the ECEF states: where ⃗ and ⃗ denote the position and velocity, respectively, in the ECEF frame.  Figure 8 shows the entire satellite covariance after 23 h in the four different frames (ECI, ECEF, RSW, and RAC). In Figure 8a and Figure 8b, the axis-to-axis correlation of a satellite is different from that of other satellites. On the other hand, all covariance ellipsoids of the RSW and RAC appear more uniform than those of ECI and ECEF. The estimated covariance in the RSW and RAC coordinate systems appear in the form of a long ellipsoid in the along-track direction. For the remainder of the study, we focused on the analysis to RSW and RAC frames.

Comparison of Orbit Correlation between RSW and RAC frames
Next, we attempted to examine the orbit correlation of each axis. Figure 9 shows the orbit correlation of each satellite and each axis corresponding to RSW and RAC at 23:00 on 15 January 2018. The blue points indicate the correlation of the RSW frame and the red points indicate that of the RAC frame. In the cases of R-S and R-A, there is a negative correlation. The along-track-cross-track and radial-track-cross-track values appear to exhibit little correlation and only in the RSW frame.

Comparison of Orbit Correlation between RSW and RAC frames
Next, we attempted to examine the orbit correlation of each axis. Figure 9 shows the orbit correlation of each satellite and each axis corresponding to RSW and RAC at 23:00 on 15 January 2018. The blue points indicate the correlation of the RSW frame and the red points indicate that of the RAC frame. In the cases of R-S and R-A, there is a negative correlation. The along-track-cross-track and radial-track-cross-track values appear to exhibit little correlation and only in the RSW frame.
In Figure 9, the along-track-cross-track correlation and radial-track-cross-track correlation of RSW are larger than those of the RAC frame. This appears to be caused primarily by the velocity change effect of the ECEF rotation. The transformation from ECI to ECEF results in a velocity change of different magnitudes depending on the location, which are A and C axes. Thus, the axis-to-axis correlations of the RSW frame differ from those of the RAC frame. In Figure 9, the along-track-cross-track correlation and radial-track-cross-track correlation of RSW are larger than those of the RAC frame. This appears to be caused primarily by the velocity change effect of the ECEF rotation. The transformation from ECI to ECEF results in a velocity change of different magnitudes depending on the location, which are A and C axes. Thus, the axis-to-axis correlations of the RSW frame differ from those of the RAC frame.

Time History of Orbit Covariance Correlation of RSW
The time history of correlation in the RSW frame is shown in Figure 10. The figure shows the time history of 29 satellites during the last 24 h for every 1-h interval. The R-S correlation is negative, while the other correlations are close to zero.

New RAC Frame Reducing R-C and A-C Correlations
The RSW characteristics are sufficient to represent the covariance. A covariance provider can represent covariance information without R-W and S-W correlation information. Actually, it is most appropriate to generate RSW coordinates by transformation from ECEF to ECI. However, the user will be subjected to a computational burden and will require the EOP. Therefore, we propose a new RAC coordinate system that does not require the EOP and has characteristics similar to the RSW. The new coordinate system is similar to the RAC frame; however, it uses a different velocity vector that negates only the rotating effect. The new coordinate system and the new velocity vector are defined as follows:

Time History of Orbit Covariance Correlation of RSW
The time history of correlation in the RSW frame is shown in Figure 10. The figure shows the time history of 29 satellites during the last 24 h for every 1-h interval. The R-S correlation is negative, while the other correlations are close to zero. In Figure 9, the along-track-cross-track correlation and radial-track-cross-track correlation of RSW are larger than those of the RAC frame. This appears to be caused primarily by the velocity change effect of the ECEF rotation. The transformation from ECI to ECEF results in a velocity change of different magnitudes depending on the location, which are A and C axes. Thus, the axis-to-axis correlations of the RSW frame differ from those of the RAC frame.

Time History of Orbit Covariance Correlation of RSW
The time history of correlation in the RSW frame is shown in Figure 10. The figure shows the time history of 29 satellites during the last 24 h for every 1-h interval. The R-S correlation is negative, while the other correlations are close to zero.

New RAC Frame Reducing R-C and A-C Correlations
The RSW characteristics are sufficient to represent the covariance. A covariance provider can represent covariance information without R-W and S-W correlation information. Actually, it is most appropriate to generate RSW coordinates by transformation from ECEF to ECI. However, the user will be subjected to a computational burden and will require the EOP. Therefore, we propose a new RAC coordinate system that does not require the EOP and has characteristics similar to the RSW. The new coordinate system is similar to the RAC frame; however, it uses a different velocity vector that negates only the rotating effect. The new coordinate system and the new velocity vector are defined as follows:

New RAC Frame Reducing R-C and A-C Correlations
The RSW characteristics are sufficient to represent the covariance. A covariance provider can represent covariance information without R-W and S-W correlation information. Actually, it is most appropriate to generate RSW coordinates by transformation from ECEF to ECI. However, the user will be subjected to a computational burden and will require the EOP. Therefore, we propose a new RAC coordinate system that does not require the EOP and has characteristics similar to the RSW. The new coordinate system is similar to the RAC frame; however, it uses a different velocity vector that negates only the rotating effect. The new coordinate system and the new velocity vector are defined as follows: where ω ⊕ denotes the angular velocity of the Earth. The error correlation of the new RAC frame is shown in Figure 11. The orbit correlation of the new RAC is similar that of RSW at 23:00 on 15 January 2018. where ω ⨁ denotes the angular velocity of the Earth. The error correlation of the new RAC frame is shown in Figure 11. The orbit correlation of the new RAC is similar that of RSW at 23:00 on 15 January 2018.

Covariance Parameterization Methods
In this section, we compare six types of covariance parameterization methods. Three are implemented by referring to an existing method or frame. Methods 1 and 2 provide the 1D RMS and 3 diagonal terms per satellite, respectively. These approaches are a reflection of how the present IGS ultra-rapid orbit generates its orbit covariance. The third method yields the diagonal components of the RAC frame, which is the frame used by IGS RTS or the PPP commercial services. Methods 2 and 3 require three parameters per satellite. In addition, we proposed three more methods using the covariance characteristics. Method 4 uses the new RAC frame. Merely by negating the coordinate rotation effect, two correlations can be eliminated. Method 5 adds the R-A correlation to Method 4. This implies that the new RAC frame neglects the R-C and A-C correlations. The last method seeks to generate complete covariance. In fact, the total covariance can be generated using six parameters: three diagonal terms and three correlation terms. Each method provides a coefficient of size to cover the actual covariance in order to take full account of the bias error caused by the orbital error.
This study evaluates each candidate method by means of the ratio of the provided covariance volume to the full covariance volume. Figure 12 shows the full and provided covariance at 24:00 on 16 January 2018. Figure 13 and Table 9 present the number of parameters required for each satellite for each method and the visualization of the full and provided volume, respectively. Table 9 summarizes the error covariance for all satellites. These results indicate that Method 1 requires only one parameter to generate the error covariance. However, it broadcasts ~ 8 times the error covariance, as indicated in Table 9. With the use of the standard deviation values for each axis, Methods 2 and 3 are able to generate 4.5 and 2.9 times the error covariance. Method 4, which uses the new RAC frame, reduces the error covariance volume to 2.2. Methods 5 and 6 (with its one additional parameter) exhibit very similar performances.

Covariance Parameterization Methods
In this section, we compare six types of covariance parameterization methods. Three are implemented by referring to an existing method or frame. Methods 1 and 2 provide the 1D RMS and 3 diagonal terms per satellite, respectively. These approaches are a reflection of how the present IGS ultra-rapid orbit generates its orbit covariance. The third method yields the diagonal components of the RAC frame, which is the frame used by IGS RTS or the PPP commercial services. Methods 2 and 3 require three parameters per satellite. In addition, we proposed three more methods using the covariance characteristics. Method 4 uses the new RAC frame. Merely by negating the coordinate rotation effect, two correlations can be eliminated. Method 5 adds the R-A correlation to Method 4. This implies that the new RAC frame neglects the R-C and A-C correlations. The last method seeks to generate complete covariance. In fact, the total covariance can be generated using six parameters: three diagonal terms and three correlation terms. Each method provides a coefficient of size to cover the actual covariance in order to take full account of the bias error caused by the orbital error.
This study evaluates each candidate method by means of the ratio of the provided covariance volume to the full covariance volume. Figure 12 shows the full and provided covariance at 24:00 on 16 January 2018. Figure 13 and Table 9 present the number of parameters required for each satellite for each method and the visualization of the full and provided volume, respectively. Table 9 summarizes the error covariance for all satellites. These results indicate that Method 1 requires only one parameter to generate the error covariance. However, it broadcasts~8 times the error covariance, as indicated in Table 9. With the use of the standard deviation values for each axis, Methods 2 and 3 are able to generate 4.5 and 2.9 times the error covariance. Method 4, which uses the new RAC frame, reduces the error covariance volume to 2.2. Methods 5 and 6 (with its one additional parameter) exhibit very similar performances.

Number of parameters per satellite
Ratio of provided covariance volume to actual covariance volume

Number of parameters per satellite
Ratio of provided covariance volume to actual covariance volume

Discussion
In this study, we analyzed real-time orbit covariance and proposed a new covariance parameterization method for low-cost user systems. Current real-time orbits provide their standard deviation without considering the correlation of each axis in the ECEF frame. Therefore, we analyzed the effect of correlation to provide a novel covariance parameterization method.
We estimated the real-time GPS orbit and covariance using DDCP observations to analyze real-time correlations. The orbit and covariance were validated using IGS final orbits. The orbit converges to the 2 cm level in the radial direction and the 7.8-cm level in terms of the 3D error after 24 h. In Figures 6  and 7, the PDF and CDF bounding plots guarantee the conservative distribution of orbit errors for each axis. The proper covariance information can be used for fault detection or user integrity.
The characteristics of the estimated covariance was analyzed over time with different frames. In Figure 7, the covariance of 29 satellites appears in the along-track direction. The covariance ellipsoid of the RSW frame appears more uniform than that of the RAC frame. The errors and covariance in the along-track direction are greater than those in the radial and cross-track directions, mainly owing to orbit dynamics. Previous studies [45] determined that the along-track error was larger than the others because it diverges in proportion to the square of the time step, whereas the errors along the other directions are only proportional to the time step. In addition, we estimated that the negative correlations observed between the radial and along-track directions are also mainly due to orbit dynamics. Positive radial error yields smaller gravity, which leads to negative along-track error [9,10]. The correlation of the cross-track and along-track errors or cross-track and radial-track errors could be neglected in the RSW frame. Although long-term characteristics [19] demonstrate that the cross-track error in the RAC frame appears unrelated to other axis errors, we determined that real-time users must consider the error correlation of each axis when using the RAC frame.
For removing the correlation of each axis, the RSW frame is the best for providing covariance. In contrast, to utilize the RSW frame, a user is required to overcome complex computation for the transformation from ECEF to ECI and additional communication for EOP parameters. Therefore, based on our analysis of covariance characteristics, we proposed a new RAC frame using the ECEF position and velocity vectors to provide covariance similar to that of the RSW frame. The new RAC coordinate uses a new velocity vector that eliminates the velocity effect due to frame rotation. Therefore, users do not have to sustain the complex computation for the transformation from ECEF to ECI in order to utilize the advantage of the RSW frame.
Finally, we evaluated six covariance provision methods that were implemented in the ECEF frame. Each method was evaluated by means of the ratio of the ellipsoid provided to the actual covariance ellipsoid, since an ellipsoid closer to the actual covariance yields more appropriate information to improve on practical applications. Method 1 used only one scalar component per satellite, and it yielded an eight-times-larger ellipsoid than the actual covariance one. Methods 2 adopted by IGS and Method 3, which utilized the standard deviations of the X-, Y-, Z-or R-, A-, C-axes, yielded nearly three or four times the full covariance. Method 4 also provided the same number of broadcasting parameters as methods 2 and 3. With the implementation of the new RAC frame, Method 4 yielded a covariance ratio of 2.2. Furthermore, Method 5 afforded a covariance ratio of 1.3 upon application of the R-A correlation. The novel proposed method is confirmed to be effective for providing covariance to users. This approach will improve covariance applications such as fault detection, integrity, and navigation performance improvement. Furthermore, applications of covariance will exhibit more continuity owing to the reduction in uncertainty. Finally, even low-cost user systems could apply the covariance information owing to the low computational burden.

Conclusions
We analyzed the characteristics of real-time orbit covariance to devise a new covariance parameterization method. For the covariance analysis, we implemented a real-time orbit determination tool. The filter in our approach utilized DDCP measurements to determine the satellite orbits and their covariance. The orbit accuracy was to the 2 cm level along the radial direction and to the 7.8-cm level in terms of 3D error relative to the IGS final orbits. In addition, we identified that the covariance conservatively reflects the error distribution.
The covariance of real-time GPS orbits exhibits a special feature in the ECI coordinate system. The error correlation between the radial and along-track errors remains negative, and the cross-track errors are uncorrelated with the other two parameters. However, the cross-track characteristics of the ECEF are variable, and they do not appear in the ECI frame or in any long-term analysis of previous study results. Therefore, we proposed a new ECEF-based local the coordinate to maintain real-time characteristics in ECI.
The new RAC frame obtained through covariance analyses is deemed suitable to provide more realistic covariance than previous approaches. The method of neglecting the correlation of each axis in ECEF yields an ellipsoid with approximately 4 times the volume of full covariance. However, the new RAC frame reduced the size of the provided ellipsoid by 2.2 times the full covariance volume. In addition, it was possible to generate an ellipsoid similar to the full covariance considering the R-A correlation. In conclusion, the proposed method provides covariance similar to the actual value with a reduced number of parameters considering the real-time covariance characteristics. We expect that the real-time covariance of the filter can be applied to navigation improvement, user integrity, and fault detection for PPP or RTK.