Abstract
The modeling of solar radiation pressure is the most important issue in precision GNSS orbit determination and is usually represented by constant and periodic terms in three orthogonal axes. Unfortunately, these parameters are generally correlated with each other due to overparameterization, and furthermore, the correlation does not remain constant throughout a long-term period. A total of 500 weeks of GPS daily solutions were estimated with the empirical CODE orbit model (ECOM) to cover various block types of satellites. The statistics of the postfit residuals were analyzed in this study, which shows the dominant annual variation of the correlations over time. There is no significant difference between eclipsing and noneclipsing satellites, and the frequency of the correlation exactly corresponds to the GPS draconitic year. Based on the residual analysis, the ECOM is the most appropriate for the Block IIR/IIR-M satellites but does not properly account for the behavior of either older Block IIA or newer IIF satellites. In addition, the daily mean residuals show a different pattern for satellite orbital planes. Therefore, the orbit model should be customized for the block types and orbital plane for better representation of multi-GNSS orbits.
1. Introduction
Since the early effort of Fliegel et al. [1] with a ROCK model, considerable progress has been made in the modeling of the solar radiation pressure (SRP) on GPS satellites. Colombo [2] introduced additional parameters to correct the resonant effects in the GPS ephemerides, which are in the form of simple empirical acceleration. The approach was adapted by Beutler et al. [3] in the Bernese GPS Software with slightly different axes to absorb the unmodeled forces in the acceleration of the satellites. This model, called the empirical CODE orbit model (ECOM), is composed of nine parameters in three orthogonal axes including the satellite-Sun vector and the spacecraft’s solar panel axis. Two periodic terms in each direction with a frequency of once per revolution are modeled under the assumption that the same errors will occur in a specific region of the orbit. The ECOM has considerably contributed to the improvement of SRP modeling, and thus, this model and its variants are extensively tested for GPS orbit determination, see Springer et al. [4], Chen and Wang [5], and Sibthorpe et al. [6], etc. Most of the International GNSS Service (IGS) analysis centers (ACs) provide the final orbit based on the adopted ECOM-type model as can be seen in Table 1.
A similar approach was made by the Jet Propulsion Laboratory (JPL) group, that is, the GPS Solar Pressure Model (GSPM) series. Following up the GSPM.II.97 model, Bar-Sever and Kuang [9] introduced an empirical model based on 4.5 years of in-flight GPS orbit data (GSPM.04), and a follow-up model for eclipsing satellites was developed [10]. The GSPM models formulated the SRP as a truncated Fourier expansion in a GPS satellite body-fixed frame, which is dependent on the Earth-satellite-Sun angle. In this model, it turns out that the estimated force parameters show a periodic characteristic in beta () angle that is the elevation angle of the Sun above the satellite orbital plane.
Despite the high performance of ECOM, a particular pattern was observed in the Satellite Laser Ranging (SLR) residuals as well as several centimeters of bias, mostly from the radial component, in IGS orbit due to Earth’s radiation pressure [11]. In addition, a box-wing satellite model estimates the rotation lag that is derived from the misalignment of the solar panels. It also adjusts the optical properties of the satellite’s surfaces and the so-called -bias, resulting in a difference at the 1-2 cm level of orbit [1, 12]. However, the partial derivatives of this model are correlated with each other (consequently the correlations in the estimated parameters) [4]; thus, it needs to constrain the parameters to the a priori values. Montenbruck et al. [13] demonstrated a new approach of SRP modeling to reduce the peaks of the radial orbit errors. It combines the simple CODE model with the cuboid-shaped body of the satellite to properly describe the acceleration, which especially improved the orbit estimation of the Galileo system.
As mentioned above, since the SRP can be decomposed into three orthogonal axes, the empirical SRP parameters are strongly correlated with each other due to overparameterization; thus, the original 9-parameter ECOM is no longer used for many years. Many ACs generally estimate five adjustable SRP terms (called “reduced” ECOM) along with the 6 initial conditions (see Table 1), that is, position and velocity at the initial epoch [14, 15]. The reduced ECOM (5-parameter) has shown a better Root Mean Squared Error (RMSE) accuracy in its orbit solution expressed in single-difference observation and the orbit overlap. However, the final model was rather empirically chosen from various combinations of SRP parameters based on a limited amount of data in the 1990’s [4, 16].
Therefore, we analyzed the GPS orbit model in terms of the correlation between SRP parameters and the residuals of orbit estimation. Specifically, almost 10 years of GPS orbit data were processed to calculate the correlations and their variation over time. The correlation was calculated from the normal matrix of the orbit estimation process. Once the orbit parameters are estimated, the residuals were analyzed to find out the contribution of each parameter in the residuals. In addition, it was also investigated to verify the characteristics of the modeling of the GPS orbit with respect to the orbital plane and the block types in the long-term series.
2. Methodology of Data Processing
2.1. Data Collection
For orbit estimation, we generated a total of 500 weeks of GPS orbits from GPS week 1421 to 1920 (DOY 091, 2007, through DOY 303, 2016). The published IGS precise orbits were used as external information to generate the orbits. Figure 1 shows the availability of satellites during the study period. Since a satellite can be replaced by a new one with the same PRN near the end of its design life, it is plotted by the satellite vehicle number (SVN). A total of 16 satellites were actively operating for the entire period, and the average lifetime of each satellite is about 6.8 years.
Figure 2 shows the block types of GPS satellites used in this study, which includes four different types, that is, Blocks IIA, IIR, IIR-M, and IIF. As can be seen in the figure, the number of IIA satellites keeps decreasing, while twelve IIR satellites are available throughout the entire period. The Block IIF satellite was first launched in mid-2010 and rapidly replaced IIA satellites. A total of 48 different satellites of almost all block types were used to analyze the parametric correlations, and the types and number of satellites are summarized in Table 2.
2.2. Solar Radiation Pressure Model
SRP modeling has been one of the most frustrating tasks in dynamic orbit estimation due to the complex properties of the satellite and the unmodeled forces acting on the satellite. Since the surface properties of the satellite and its attitude about the reference axes are not known with sufficient accuracy, the empirical parameters should be solved in the orbit determination process [17]. Although many SRP models were proposed and applied by ACs, the ECOM including its variants was the most successful model in GNSS orbit determination [3]. The model is basically composed of three orthogonal directions with constant and once per revolution terms in each direction [4, 16]. The principal direction is along the vector of the satellite-Sun direction, the second axis points to one of the solar panel axes of the satellite, and the last axis assures the orthogonality: where is the a priori acceleration from a model, is the argument of latitude of the satellite, and and represent the unit vector of the satellite-Sun direction and the vector along the spacecraft solar panel axis, respectively. It assumes a nominal, yaw-steering satellite attitude, and lastly, completes the right-handed system defined by . The a priori model can be obtained from the previous model [16], although in the case of the CODE analysis center, no a priori model has been used since July 2013 [18] due to the lack of reliable coefficients for new GPS and all GLONASS satellites.
In practice, the accurate control of satellite attitude in space is difficult; thus, the empirical parameters should be estimated. The original ECOM is known to be well suited for the near cubic-shaped GPS satellites but not for the elongated cylinder-shaped GLONASS-M satellites. It uses constant terms along with the once per revolution parameters in each direction, which is represented by sine and cosine terms as a function of the argument of the latitude on the satellite orbital plane. However, the absence of periodic terms in is regarded as a major deficit of the reduced ECOM. Since January 4, 2015, ECOM has been extended to include different frequencies of the argument for periodic terms in - and -axes [15] with the assumption of perfect attitude control of the satellite, as given by the following:
Contrary to the original ECOM, the new extended ECOM provides information based on twice per revolution as well as fourfold per revolution terms in some direction. In addition, the even-order terms are available only in the component, while the odd-order periodic terms are modeled on the -axis. The arguments of the periodic terms, , are slightly different from the original model, that is, the argument of latitude for the satellite with respect to that of the Sun. It provides even better intuitive formulation of the estimated parameters although the differences are negligible for short one-day arcs [15]. One thing to be mentioned is that the acceleration due to SRP should be turned off when the satellite is in eclipse or scaled down according to the fractional area of the Sun as seen from the satellite. The new extended ECOM was applied by CODE since early 2015, and therefore, the original ECOM was used in this study throughout the entire period to ensure the consistency of the SRP model.
2.3. Orbit Integration
The forces acting on the satellite are generally a function of the position and velocity of the satellite. The calculated accelerations are integrated to generate the states of the satellites at later time epochs. For the reconstruction of the orbit, we developed a bidirectional, multistep numerical integrator for the acceleration. A total of about 10 years of GPS orbit data were integrated for the analysis. Table 3 shows the summary of the orbit integration strategy. Most of the models refer to the International Earth Rotation and Reference Systems Service (IERS) Conventions 2010 [19].
2.4. Orbit Validation
It is very important to evaluate the orbit solutions because the deficiencies in SRP modeling can be assessed through this process. There are several ways to check the orbit quality, which are most commonly used in the orbit validation: (1) the postfit residuals for internal consistency, (2) the misclosures at midnight epochs, (3) the satellite laser ranging as an independent validation of the orbit, (4) the estimates of the geocenter coordinates, (5) the Earth Rotation Parameter (ERP) differences with respect to the IERS solutions, and (6) the scale parameters of Helmert transformations for station coordinates [15, 21, 22].
The orientation of the orbit can be entirely determined by the IGS orbit as pseudoobservations. However, it is necessary to verify the quality of the orbit solution with completely independent external information, that is, SLR data. In addition, the addition of SLR observations in GNSS orbit determination does not improve accuracy due to the low availability of the SLR data [23]. It is reported that there are no scale issues in the GNSS and SLR technique solution (consistent at the 1 mm level) model [24], although there are some biases due to the SLR receiver systems. Nevertheless, only two GPS satellites were equipped with SLR reflectors during most of the period. Therefore, we used the postfit residuals to validate the performance of the orbit modeling in this study, although the residuals depend on a priori values to some extent.
3. Experiments
3.1. Cross-Correlation of ECOM Parameters
For the analysis of the SRP modeling performance, we estimated all 9 parameters (3 components of constant/cosine/sine terms in each axis of the frame, see equation (1)) from the original ECOM. Figure 3 shows the indices of the correlation between SRP parameters. Since the autocorrelation should be 1, the diagonal terms are omitted intentionally. Therefore, a total of 36 correlations were analyzed to check the variation of the correlation in time series.
Figure 4 shows the exemplary daily correlation coefficients averaged for all satellites in absolute sense on DOY 195, 2009, and the autocorrelation was excluded for convenience. The correlation should be, by definition, either positive or negative. However, since we need to analyze the magnitude of the correlation for each parameter combination, and the correlation may be changed day after day, only the average correlations in absolute sense were considered in this study. One thing to be noticed is that the averaged magnitude becomes almost twice with the absolute values, which can be seen in Figures 4 and 5.
As can be seen in the figure, the parameters in the - and -axes are highly correlated with each other (see and ), and the direct Sun-satellite component also correlated as well. These patterns are distinctly consistent throughout the entire period, although there are some differences depending on the satellite block types and individual satellites.
Figure 5 represents the correlation coefficients for all combinations of SRP parameters for eclipsing and noneclipsing conditions, where the indices refer to Figure 3. As discussed above, the magnitude of the correlation in absolute sense is almost twice the magnitude of actual correlation (either positive or negative). It shows that there is no significant difference due to eclipsing states for most of the parameters. This is partially because the ACs use essentially the same model to obtain the orbits; thus, the eclipsing satellites perform equally well.
However, the maximum correlations can be observed for different parameter sets; for example, the eclipsing satellites have a maximum correlation for and (indices 10 and 22, respectively). The sine and cosine terms in the direct satellite-Sun direction are correlated with each other, and the nominal yaw attitude does not seem to work properly while in eclipsing states. On the other hand, the maximum correlation happens to and (also indices 28 and 34, respectively) for the noneclipsing satellites, which is attributed to the fact that the - and -axes are orthogonal to each other and the sine/cosine terms are shifted 90 degrees in phase.
Table 4 shows the number of satellites experiencing eclipses due to the shadow of the Earth, along with the duration in minutes. Overall, about 7 satellites are in eclipse conditions each day, reaching a maximum of 12 per day, and the average duration is over 40 minutes.
Figure 6 shows the different types of correlation coefficients between SRP parameters. As can be seen in the figure, the correlation between the - and -axis parameters is almost constant with only slightly different magnitudes (Figures 6(a) and 6(c)), although there are fluctuations with a combination of several frequencies of the signal. This indicates that there is a constant correlation between two axes, and these correlations do not show a temporal variation in the long-term series. It is interesting that the -axis shows a temporal periodic variation of almost annual motion with respect to both the - and -axes (Figures 6(b) and 6(d)) but slightly different characteristics in their patterns. The - and -axes are highly correlated with each other with an apparent frequency of the signal, while the sine and cosine terms in the - and -axes show a temporal variation with a mixture of different signals.
(a)
(b)
(c)
(d)
From Figure 6(d), the annual variation in the correlation coefficient between and is clear for the entire period. To find out the frequency of the signals, this correlation was transformed into the frequency domain, resulting in the dominant once per revolution frequency (see Figure 7). Other than this frequency, three revolutions per year are also significant but with a much smaller amplitude. It is known that a perturbation in the GPS satellite orbits occurs at harmonics of the GPS draconitic year due to the relative position of each satellite, Earth, and Sun [11, 14]. The most dominant frequency in Figure 7, that is, about 1.04 cycle per year, corresponds to the GPS draconitic year as pointed out by Ray et al. [25]. Therefore, the residual analysis of the orbit model also supports the harmonic signals in the time series of IGS orbits.
3.2. Residual Analysis
The orbit representation of this study refers to the general Gauss-Markov model. The IGS final orbits were used as external information; thus, the difference with the obit propagated from the initial conditions becomes an observation. The design matrix is the partial derivatives of the observation with respect to the unknown parameters, which cannot be calculated directly due to the implicit relationship, as given in where and represent the actual and the computed measurements, respectively; and denote position and velocity vectors; is the dynamic parameter of SRP; and is the random measurement error.
Thus, the design matrix of the Gauss-Markov model can be represented by the variational partials which are simultaneously integrated with the acceleration by the newly developed numerical integrator. The overall RMSE of the postfit residuals is calculated based on the following equation [7]: where is the residuals after the least squares adjustment, is the number of satellites to be integrated, and represents the total number of observations in the calculation. Since the estimated orbit uses similar dynamic models to the IGS final orbit, it can be expected to generate small residuals, although not all ACs use the same model and the stochastic pulses.
Once the propagated orbit was fitted to the published one, the RMSE of the residuals can be calculated for each day. Figure 8 represents the mean and standard deviation of the residuals each day, plotted separately, for the entire period of this study. The mean values are almost constantly zero but apparently fluctuate with an amplitude of about ±1 mm at a frequency of once per year. The first Block IIF satellite was launched in 2010 (denoted as (1)), and the number of IIF satellites has overtaken that of IIA satellites as of mid-2014 (denoted as (2)) as inferred from Figure 2.
On the other hand, the standard deviation of the residuals shows a different behavior from the mean value over time. The entire period can be divided into three sections by the vertical lines (1) and (2). In the first part, before mid-2010, the residuals decrease almost linearly. However, the RMSE became considerably stable since Block IIA satellites started being replaced by Block IIR-M satellites in 2009. This trend continues until mid-2014 when Block IIF satellites began to replace more IIA satellites. Since the year 2014, in which Block IIA satellites make up only a small part of the satellite constellation, the RMSE of the residuals has almost bounced back up to the level seen at the beginning of the test period. Therefore, since Block IIR satellites are predominant during the second section (between (1) and (2)), it can be concluded that the ECOM is best fitted to Block IIR/IIR-M satellites.
3.3. Effect of Eclipses and the Orbital Plane
The satellites experience an eclipse when sunlight is blocked by the shadow of the Earth or Moon, depending on the geometry with respect to the satellite. Eclipses caused by the Earth mainly occur when the beta angle is less than about ±15 degrees, while lunar eclipses happen regardless of the satellite orbital plane. Figure 9 shows the RMSE of the residuals with respect to the beta angle for 3,500 days of all satellite orbit solutions. Contrary to expectations, the RMSE is slightly larger around the beta angles of ±40 degrees but not significantly so. For a clearer interpretation of the values, moving average RMSE values are plotted every 5 degrees (in terms of beta angles) with large yellow circles. The result agrees with that of Figure 5; that is, there are no significant degradations in the orbit solution caused by eclipses.
Since the satellites in the same orbital plane show a similar behavior, we checked the residuals by the orbital plane (see Figure 10). The residuals of each satellite were grouped into the orbital plane and averaged each processing day. The mean values by each plane are plotted with an offset for the entire period as denoted in the figure. Similar to Figure 8, the residuals show a large fluctuation on each side but have different patterns for the orbital plane. In particular, the residuals of the satellites in plane D are considerably stable since mid-2011, while those in plane B keep increasing for the period in discussion. Therefore, it may be feasible to model the SRP parameters considering the characteristics of the orbital planes.
4. Summary and Conclusions
SRP modeling has been the most important issue for several decades in satellite orbit determination. Many SRP models were proposed, of which the ECOM was the most successful model. However, the correlation between the SRP parameters is still a difficult problem to resolve. We calculated 500 weeks of GPS orbit solutions (year 2007 to 2016) to validate the SRP models based on the correlation of the parameters and the residuals. Most block types of GPS satellites were analyzed during this period including 16 full-time IIR satellites.
Since the orbit solutions by IGS ACs have based on the ECOM-type SRP parameterization until very recently, we also adopted this system throughout the study period to avoid any constraints on specific parameters. The Gauss-Markov model was applied to estimate the parameters, and the variational partials were integrated for the design matrix along with acceleration. Although SLR observation is often used as independent external information for orbit validation, the residuals of the orbit solution were analyzed in this study. This is because the main purpose of this study is to analyze the correlation between the SRP parameters, and the amount of SLR observations is not sufficient as well. All correlations between the original ECOM parameters, excluding the autocorrelation, were analyzed in time series.
The components in the direct Sun-satellite direction are correlated with most of periodic terms, while the - and -axes are highly correlated with each other. The eclipsing condition does not impair the orbit residuals for most of the correlations. However, the maximum correlation was observed for different combinations of SRP parameters. That is, the - and -axes are highly correlated for noneclipsing satellites, while they are related to the direct Sun-satellite direction for the eclipsing satellites. This suggests that there may be some problems in nominal yaw attitude control, and it is necessary to include yaw modeling during eclipses.
The correlation between the - and -axis parameters is nearly constant in time but shows a small fluctuation of several different frequencies. On the other hand, the -axis is highly correlated with the - and -axes with an apparent annual variation which precisely corresponds to the GPS draconitic year. Therefore, it can be mentioned that the residual analysis of the orbit model also supports the harmonic signals in the time series of IGS orbits. The time series behavior of the residuals shows the best performance during the midsection (2010 to 2014) at which Block IIR/IIR-M satellites are dominant in number. Therefore, the ECOM seems to work better for Block IIR/IIR-M satellites, while other parameterization may be necessary for Block IIA and/or IIF satellites. In addition, there is no significant difference by the eclipsing condition as well as the geometry of the orbital plane with respect to the Sun. However, the daily mean residuals show a different pattern for each orbital plane, which should be considered for the customized parameterization of the SRP.
It was reported that the addition of twice per revolution terms in the -axis improves the quality of the orbit. New extended ECOM reduces the peaks of the draconitic year harmonics in GNSS-derived Earth Orientation Parameters (EOPs) and the day boundary misclosures [15]. Some research argues that multifrequency per revolution can help the improvement of the quality of the orbit solutions. Therefore, it is necessary to understand the correlation between the parameters of higher frequencies, resulting in customized SRP models for each block types in further study.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2018R1D1A1B07048475).