Introduction

The US global positioning system (GPS) has been providing positioning, navigation and timing (PNT) services for more than 20 years (Johnston et al. 2017). Since the launch of the first Block I vehicle, GPS satellites have been developed from generation to generation (Block I, II, IIA, IIRA/B/M, IIF, III). The last Block I satellite was taken out of service in November 1995, while the final Block II satellite was decommissioned in March 2007. The first new Block III satellite was launched in December 2018 and was unusable from July to October 2019. The second and the third GPS III satellites were launched in August 2019 and June 2020, respectively (https://www.navcen.uscg.gov). Currently, the data period of Block III satellites is too short to analyze this study. Thus, this paper focuses mainly on GPS Block IIA, IIR and IIF satellites.

The current IGS final GPS orbit products have a 1D mean RMS value over the three geocentric components of about 2.5 cm based on comparisons with independent laser ranging results and discontinuities between consecutive days (http://www.igs.org). One of the most important errors in current GPS orbit models is due to the effect of solar radiation pressure. As Springer et al. (2019) presented, considerations of SRP modeling for GPS satellites in individual analysis centers are not identical and performances of each model for each Block type are different as well. Consequently, it is hard to conclude which SRP model is better by making comparisons of orbit solutions between different analysis centers.

SRP can be modeled both analytically and empirically. Milani et al. (1987) formulate the physical interaction between SRP and satellite metadata, such as satellite attitude, mass, dimensions and optical properties. The operation law of the current GPS satellite attitude is clearly known except for that of Block IIA satellites inside the eclipsing season. Metadata of GPS Block II/IIA and IIR satellites are published in great detail by Fliegel et al. (1992) and Fliegel and Gallini (1996). The shape, area and optical properties of each element on each satellite surface are given individually. For the sake of simple usage, Rodriguez-Solano (2014) computes the area of each satellite surface. Optical properties of each element are averaged into surface-specific properties according to the proportion of their individual areas. In addition, Rodriguez-Solano (2014) gives dimensions and optical properties of GPS Block IIF satellites based on some unpublished documents. For the latest GPS Block III satellites, metadata is published by (Steigenberger et al. 2020). Since metadata of satellites cannot be known without any error and some other perturbations contribute as well, it is nearly impossible to model SRP perfectly by using an analytical model.

In the absence of a precise analytical model, empirical ROCK4 and ROCK42 models are first developed by the satellite manufacturer Rockwell International for Block I and II satellites (Fliegel et al. 1985; Fliegel and Gallini 1989). Both ROCK models can be represented as short Fourier series as a function of the elongation angle (Earth–Sun–satellite) in X and Z directions (satellite body-fixed frame) while neglecting the Y bias. The standard ROCK models (version “S”) are later enhanced considering the thermal model and are named as version “T.” As an extension of the ROCK models, Bar-Sever and Kuang (2004) from JPL (Jet Propulsion Laboratory) developed a new GSPM (GPS solar radiation pressure model) model for GPS satellites. This model is later extended for GPS satellites inside eclipse seasons as well (Bar-Sever and Kuang 2005). Compared with the ROCK models, the GSPM model employs more Fourier terms as a function of elongation in the X direction. In addition, the Y bias is modeled as well. All the GSPM parameters are derived by fitting precise JPL GPS satellite orbits with an arc length of 10 days over several years. As presented by Sakumura et al. (2017), the GSPM model is updated and extended to Block IIF satellites by making use of the latest precise GPS orbit products. However, detailed parameters are not given.

Both ROCK and GSPM models are developed based on the analysis of satellite accelerations that are caused by SRP or radiation-related perturbations. Knowledge of satellite metadata is a good precondition. Due to the lack of such information, the ECOM model is developed in the early 1990s (Beutler et al. 1994). In the nominal yaw steering mode, the ECOM parameters are defined in a Sun-oriented DYB frame, with D pointing toward the Sun, Y along the solar panel axis and B completing the right-handed system. The classical ECOM model consists of three constant parameters (D0, Y0 and B0) and six first-order Fourier coefficients (DC, DS, YC, YS, BC, BS) in the D, Y and B directions. Then, instead of using all the nine parameters, Springer et al. (1999a) from CODE concludes that the periodical terms in D and Y directions are not necessary for GPS satellites and only the reduced 5-parameter ECOM model (D0, Y0, B0, BC, BS) is routinely used. As an extension of the ECOM model, a new ECOM2 model was developed by CODE, initially serving GLONASS satellites (Arnold et al. 2015). They find that the spurious signals in the GLONASS-derived geocenter z-coordinate are significantly reduced compared to that by using the ECOM model. Prange et al. (2017) introduce both ECOM and ECOM2 models into GNSS satellite orbit determination. The ECOM2 model performs better than the ECOM model for Galileo and QZSS (Quasi Zenith Satellite System) satellites. The reason is that both Galileo and QZSS satellites are elongated in shape. SRP accelerations change periodically due to the illumination of each satellite surface. In this case, the higher-order Fourier terms in the ECOM2 model can compensate for such periodical accelerations fairly well. However, more parameters in D direction mean potentially stronger correlations between parameters. This can weaken the solutions of other estimated parameters, for instance, the length of day (LOD).

To reduce the number of ECOM parameters or to constrain some of the ECOM parameters, an a priori model is preferred. In principle, both ROCK and GSPM models can be taken as a priori model and jointly used with ECOM models. However, ROCK models are only available for Block I and II satellites and the latest GSPM parameters are not yet published. Springer et al. (1999b) generalize an a priori model as a function of \(\beta\) (Sun elevation above the orbital plane) angle by fitting CODE precise orbits over long time periods. The model is later recomputed using CODE final orbits from the year 2000 to 2006 (Dach et al. 2009). However, this a priori model is no longer used in CODE's routine processing. Another classical a priori model is the box-wing model. As proven by Bury et al. (2020), Li et al. (2019), Montenbruck et al. (2017, 2015b), Steigenberger et al. (2015) and Zhao et al. (2018), an a priori box-wing model improves satellite orbits significantly for a satellite with an elongated shape. Therefore, the main purpose of this research is to set up an enhanced box-wing model for GPS satellites.

In the following sections, we first describe the current GPS satellite attitude and the impact on modeling SRP. Then, we show the impact of all the discussed physical effects on SRP. Third, we adjust all the box-wing parameters, including additional physical parameters, such as yaw bias, radiator emission and thermal radiation (Duan et al. 2020, 2019b; Rodriguez-Solano et al. 2012a). In addition, the GSPM parameters are determined as well by using the same GPS measurements. Fourth, both the enhanced box-wing model (using all the adjusted box-wing parameters and physical effects) and the GSPM model are taken as the a priori model and are jointly used with ECOM models in GPS satellite orbit determination. Finally, we summarize conclusions and outlooks.

GPS satellite attitude

The box-wing model depends on the satellite attitude. Also, the consideration of physical effects on satellite surfaces, in particular the radiator emission, requires accurate satellite attitude both outside and inside eclipse seasons. GPS satellites take a nominal yaw steering attitude in non-eclipsing seasons. Two requirements determine the nominal attitude: first, that the transmitting antenna points toward the geocenter, and second, that the solar panel axis is perpendicular to the Sun–satellite direction. To satisfy these two conditions, the satellite body has to yaw constantly. This task is directed by the onboard attitude control system (ACS) based on the data from solar sensors. As long as the Sun is visible, data from solar sensors could be a representation of the yaw error (Bar‐Sever et al. 1996). Then, the ACS decides the correct yaw direction to minimize this observed yaw error.

The eclipse season is the period when a satellite crosses the Earth's shadow once per revolution. For Block IIA satellites during the shadow crossing, solar sensors provide no data and the ACS might be driven in an open-loop mode by the system's noise. In order to avoid this, a small yaw bias of + 0.5 deg is artificially added for Block IIA satellites (Bar-Sever 1996). As described by Bar‐Sever et al. (1996), upon shadow entry, a Block IIA satellite either maintains or reverses its natural yaw direction with its maximum yaw rate depending on the sign of the yaw bias and the \(\beta\) angle. Upon shadow exit, the Block IIA satellite needs about 30 min to recover to the nominal orientation. However, this post-shadow recovery is hard to be correctly modeled since the optimal yaw direction is unclear. The shadow crossing behaviors of Block IIR and IIF satellites are different from that of Block IIA satellites. Block IIR satellites can maintain nominal yaw orientation during a shadow crossing by using the constant hardware yaw rate of 0.2 deg/sec (Kouba 2009). Block IIF satellites compute the yaw rate at a shadow entry, ensuring that the nominal yaw rate is reached at the shadow exit. This means no yaw angle recovery is needed at the shadow exit for Block IIF satellites (described by the IGS “eclips.f” Fortran routine, developed by J. Kouba).

Noon and midnight turn maneuvers happen when the satellite–Earth and satellite–Sun vectors are close to collinear. When the nominal yaw turn rate exceeds the respective maximum hardware rates, the GPS satellites start turning with the maximum hardware rate until the yaw angle catches up with the nominal yaw rate. This may take up to 30 min for Block IIA and IIF satellites (maximum yaw rate of 0.10–0.13 and 0.11 deg/sec for about 180 deg turn) and 15 min for Block IIR satellites (maximum yaw rate of 0.20 deg/sec for about 180 deg turn). For Block IIF satellites with small negative \(\beta\) angle (> − 0.9 deg), Dilssner (2010) observes that during noon turns, Block IIF satellites rotate in a negative direction. The reason is later confirmed by a negative yaw bias of − 0.5 deg.

In this contribution, we use the IGS “eclips.f” Fortran routine to represent the GPS satellite attitude. The satellite body frame orientation for Block IIR satellites is modified to the IGS definition (Montenbruck et al. 2015a). The post-shadow recovery periods (about 30 min) for Block IIA satellites are not properly modeled inside the routine as the exact behavior of the satellite cannot be predicted (Kouba 2009). The yaw direction (either normal or reversed rotation) selected by the satellite to optimally recover to the nominal orientation is uncertain. Since GPS Block IIA satellites have a significant antenna phase center offset in the X direction the impact of the improper attitude modeling on phase and pseudorange observations can be as large as 10 cm. The phase wind-up error also depends on the satellite attitude; the impact on undifferenced phase measurements can also reach up to 10 cm (Kouba 2009). Therefore, Block IIA satellites measurements during the post-shadow recovery periods are suggested to be excluded for precise applications.

However, the impact of the SRP cannot be avoided if an a priori box-wing model is considered. Also, we conclude in Sect. 3 and Sect. 4 that Block IIA satellites have constant radiator acceleration in satellite body-fixed − X (− 0.9 \(nm/s^{2}\)) and Y (− 0.3 \(nm/s^{2}\)) directions. To investigate the impact, we compute acceleration differences caused by SRP and radiator emission between the normal yaw direction and the reversed yaw direction for one Block IIA satellite during a post-shadow recovery period. The SRP is assumed to be fully described by a box-wing model. As shown in Fig. 1, acceleration differences in the D direction can be larger than 2.5 \(nm/s^{2}\), while in the Y direction, acceleration differences are minor. Figure 2 shows the corresponding orbit differences by fitting a 3-day-arc of CODE precise orbits, as well as differences of the orbit prediction of 1 day. No empirical parameters are estimated in the top figure, while the ECOM model is employed in the bottom figure. 3D RMS of orbit differences of the fitted 3-day-arc is about 3.5 cm for the result without any empirical parameters, while it is about 0.5 cm for the result with the ECOM model. The impact on the predicted orbits is larger, with an RMS value of 11.1 and 2.0 cm, respectively, for the result without and with the ECOM model. We can conclude that the uncertainty of the post-shadow recovery rotation direction of Block IIA satellites has a minor dynamic impact on the satellite orbit.

Fig. 1
figure 1

Acceleration differences of SRP and radiator emission in the ECOM DYB frame by using the normal and the reversed yaw direction during a post-shadow recovery period for one Block IIA satellite (\(\beta\) = 6 deg, post-shadow recovery starts from epoch 0 and lasts about 30 min)

Full size image
Fig. 2
figure 2

Impact of post-shadow recovery periods on satellite orbits for one Block IIA satellite (\(\beta\) = 6 deg). The top figure shows the result without any empirical parameters, while the bottom figure shows the result with the 5-parameter ECOM model. The solid black vertical line denotes the starting epoch of each post-shadow recovery period, while the green vertical line represents the starting time of the 1-day-arc orbit prediction

Full size image
Fig. 3
figure 3

Simulated impact of yaw bias of + 0.5 deg on SRP accelerations for one Block IIA satellite in the ECOM DYB frame, \(\beta\) denotes the Sun elevation, \(\Delta u\) represents the argument of latitude of the satellite with respect to the Sun

Full size image

Physical effects

Apart from regular SRP, a number of other physical effects contribute as well. Yaw bias can be caused by a solar sensor bias and is a function of the elongation angle. X- and Y-axes of the satellite body are rotated by the yaw bias around the Z-axis with respect to the original satellite body-fixed frame. Therefore, SRP effect of satellite surfaces may change accordingly. As described by Milani et al. (1987), SRP acceleration equations for solar panels and satellite body surfaces are different since multilayer insulation (MLI) blankets usually cover the satellite body. The energy absorbed by the satellite surfaces is assumed to be instantaneously reradiated back into space, according to Lambert’s law (Milani et al. 1987). To combine these two equations, we introduce a thermal reradiation factor \(kappa\), as follows:

$$ {\mathbf{acc}} = - \frac{A}{M}\frac{{S_{0} }}{c}\cos \theta \left[ {\left( {\alpha + \beta } \right){\mathbf{e}}_{D} + \frac{2}{3}\left( {\delta + \kappa \alpha } \right){\mathbf{e}}_{N} + 2\rho \cos \theta {\mathbf{e}}_{N} } \right] $$
(1)

where \(A\) denotes the surface area, \(M\) the total mass of the satellite, \(S_{0}\) the solar flux, \(c\) the vacuum velocity of light, \(\kappa\) the thermal reradiation factor (0 for solar panels and 1 for satellite body surfaces in this contribution), \(\alpha\), \(\delta\), \(\rho\) the fractions of absorbed, diffusely scattered and specularly reflected photons. Furthermore, \({\mathbf{e}}_{D}\) denotes the Sun direction, \({\mathbf{e}}_{N}\) the surface normal vector and \(\theta\) the angle between both vectors.

The yaw bias impact on SRP can be easily computed by using (1) considering a biased \({\mathbf{e}}_{N}\). We would like to mention that for the solar panel only the \(\delta\) and \(\rho\) parameters contribute to the yaw bias impact because \(\alpha\) does not cause an acceleration in the direction of \({\mathbf{e}}_{N}\). As solar panels are designed to absorb as much energy as possible, the corresponding \(\delta\) and \(\rho\) are usually small. Figure 3 shows a simulation of one GPS Block IIA satellite with a yaw bias of + 0.5 deg. SRP accelerations are described by the same box-wing model as that in Fig. 1. The displayed accelerations are the impact of the yaw bias on the SRP in the ECOM DYB frame. It is noted that the largest impact appears in the Y direction; the total averaged amount over the whole period is about -0.15 \(nm/s^{2}\). As discussed by Fliegel et al. (1992), the constant empirical estimates of Block II/IIA satellites in the Y direction (Y bias) are caused by two main reasons. First, the yaw bias of about + 0.5 to + 1.0 deg causes perturbations in the Y direction. Second, the heat generated in the satellite body is radiated from louvers on the Y side of the satellite. Typically, the Y bias is considered to be dominated by the yaw bias. To investigate this, we compute Y bias (Y0) estimates of one year for all the GPS satellites by using the 5-parameter ECOM model without any stochastic pulses. Figure 4 shows the daily Y0 estimates of three GPS satellites. Block IIA and IIR satellites show clearly larger Y biases than that of the Block IIF satellite. Considering the impact of the yaw bias of + 0.5 deg that we have computed, it seems likely that the radiator emission mainly causes Y0 estimates of Block IIA and IIR satellites. This radiator effect in Y direction for the Block IIF satellite is, however, either balanced on both sides or is very small.

Fig. 4
figure 4

Daily 5-parameter ECOM Y0 estimates of one year for one IIF, one IIR and one IIA satellite

Full size image

It is also noted that all the three satellites show clearly anomalous Y0 estimates inside the eclipsing season. This indicates that the radiator emission exists when a satellite crosses shadows. Since satellite attitude inside the eclipse is different from the nominal yaw, a potential radiator effect in the − X direction can be wrongly absorbed by the ECOM Y0 parameter. In fact, this provides a new approach to confirm satellite attitude inside eclipses. For this, we set up two additional constant radiator accelerations in the satellite body-fixed − X and Y directions inside shadows for each satellite and jointly estimate them together with the 5-parameter ECOM parameters in GPS satellite orbit determination (one-day-arc solution). Figure 5 shows daily radiator estimates and formal uncertainties of one Block IIF satellite in the satellite body-fixed Y direction using both the nominal attitude and the attitude described in “eclips.f”. At the edges of eclipse seasons with \(\beta\) angles close to 14.5 deg, accelerations are not well determined since the shadow periods are short and the number of observations is low. Estimates using the nominal attitude show clear anomalies when the \(\beta\) angle is close to zero and the mean formal uncertainty is also about two times larger than that of the “eclips.f” attitude. The “eclips.f” attitude provides constant estimates every day except for the day with a starting \(\beta\) angle of − 0.93 deg. The reason could be that the “eclips.f” attitude considers a negative yaw bias of − 0.5 deg when the \(\beta\) angle is greater than or equal to − 0.9 deg. However, in the real case the true threshold might be slightly more or less than − 0.9 deg. Figure 6 shows similar radiator estimates for one Block IIA and one Block IIR satellite. Values regarding the Block IIA satellite are shifted by a constant value of − 20 \(nm/s^{2}\). Differences between the nominal attitude and the “eclips.f” attitude for the Block IIR satellite are very small because Block IIR satellites can maintain nominal attitude inside eclipse seasons except for the midnight and noon turns. The “eclips.f” attitude describes these two maneuver turns correctly, and we observe more constant estimates when the \(\beta\) angle is close to zero. The Block IIA estimates perform differently than that of the Block IIR and the IIF satellite. Daily radiator estimates by using both types of attitude show clear variation. The reason could be that the employed satellite yaw rates might not be accurate enough.

Fig. 5
figure 5

Daily radiator estimates and formal uncertainties of one Block IIF satellite in the satellite body-fixed Y direction using both the nominal and the “eclips.f” attitudes

Full size image
Fig. 6
figure 6

Daily radiator estimates and formal uncertainties of one Block IIA and one Block IIR satellite in the satellite body-fixed Y direction using both the nominal and the “eclips.f” attitudes. Values regarding the Block IIA satellite are shifted by a constant value of − 20 \(nm/s^{2}\)

Full size image

Another potential physical effect is outgassing. As described by Fliegel and Gallini (1996), there are four candidate sources of gas: the solar panels, the apogee engine, the SV body and the MLI material. Soon after arriving in orbit, the SV is affected by the emission of the gas that is incorporated inside these materials. In order to show the effect, we compute the D0 parameter of one Block IIF satellite that was launched in June of the year 2014. The 5-parameter ECOM model is employed without any a priori model in the satellite orbit determination. Figure 7 shows the daily D0 estimates of about 6 years since the data have been officially available. We find that the outgassing effect lasts about half a year after the launch of the satellite. The largest effect reaches up to 2 \(nm/s^{2}\). Therefore, this effect cannot be ignored in the following optical parameter or the GSPM parameter adjustment.

Fig. 7
figure 7

Daily 5-parameter ECOM D0 estimates of one Block IIF satellite

Full size image

Adjustment of enhanced box-wing and GSPM parameters

In this section, we determine box-wing and GSPM parameters, respectively, based on the same GPS code and phase measurements covering 6 years. We use about 150 globally distributed tracking stations. The time span starts from doy (day of year) 1 2014 to doy 365 2019. The outgassing and maneuver time periods are excluded in the box-wing parameter adjustment. Measurements inside eclipse seasons are additionally excluded in the GSPM parameter adjustment to avoid a negative impact on the estimated parameters. We estimate station coordinates and tropospheric delays as a first step using PPP (precise point positioning) approach with CODE final GPS products. Then, station-related parameters are fixed in the box-wing and the GSPM parameter adjustment. The orientation of the GPS satellite body frame in this study follows the IGS definition.

Enhanced box-wing model

Since the + X surface of Block IIA and IIR satellites consists of both flat and cylindrical parts, we introduce a shape factor to combine SRP equations for both types of surface, as expressed in

$$ {\mathbf{acc}} = - \frac{A}{M}\frac{{S_{0} }}{c}\cos \theta \left[ {\left( {\alpha + \beta } \right){\mathbf{e}}_{D} + \left( {\frac{\pi }{6}s + \frac{2}{3}\left( {1 - s} \right)} \right)\left( {\delta + \kappa \alpha } \right){\mathbf{e}}_{N} + \left( {\frac{4}{3}s + 2\left( {1 - s} \right)} \right)\rho \cos \theta {\mathbf{e}}_{N} } \right] $$
(2)

where s denotes the shape factor, indicating the weighting between a cylindrical surface (1) and a flat surface (0). Optical parameters \(\alpha + \delta\) and \(\rho\) are adjusted in + X and \(\pm\) Z surfaces for each satellite (Duan et al. 2019b). The solar panel \(\delta\) parameter is estimated to compensate for the thermal radiation performed on solar panels (Duan et al. 2020). A rotation lag of solar panels is estimated for each satellite (Rodriguez-Solano et al. 2012a). Furthermore, we assume a constant solar sensor bias (causing the yaw bias) for each satellite in the adjustment and compute its numerical partial derivative. As the yaw bias effect is considered, the remaining Y bias is mainly caused by the radiator effect. So, we estimate a constant acceleration in the Y direction both inside and outside the eclipse. To compute the radiator effect from the − X surface, we have two iterations. First, we estimate additionally a constant acceleration inside shadows in the − X surface for each satellite. Then, we apply this estimated acceleration constantly both inside and outside eclipses and determine all the other parameters similarly as in the first iteration. We find that in the second iteration, the acceleration inside shadows in the − X surface is very close to zero. This is a way to reduce the correlation between the − X acceleration and optical parameters of the + X surface. Therefore, we have six optical parameters of satellite body surfaces, one optical parameter of the solar panel, one rotation lag of the solar panel, two radiator parameters and one solar sensor bias, which is a total of 11 parameters for each satellite.

The processing settings are shown in Table 1. Earth rotation parameters (ERPs) are fixed to the IERS Bulletin A values. Earth albedo (Rodriguez-Solano et al. 2012b) and antenna thrust (Steigenberger et al. 2018) are considered. Undifferenced phase ambiguities are fixed to integer values using a modified Bernese software (Dach et al. 2015). Dimensions and total mass of a satellite are fixed to the published values. Optical parameters from Rodriguez-Solano (2014) are taken as a priori values in the least-squares adjustment. Daily normal equations of 6 years are stacked after pre-eliminating all but the 11 box-wing parameters. For the satellite-specific estimates, we compute the mean value of each Block type.

Table 1 Settings of box-wing parameter adjustment
Full size table

Tables 2, 3 and 4 show all the adjusted box-wing parameters. The total radiator emission power is about 340 W, 290 W and 380 W for Block IIA, IIR and IIF satellites, respectively, considering the estimated radiator accelerations. All the formal uncertainties are reasonably small. Unphysical negative optical estimates are observed for some surfaces. The summed values of these three optical estimates of the same surface are not equal to the physical condition 1. One of the reasons is that parameters \(\alpha + \delta\) and \(\rho\) of the same surface are correlated. Systematic errors might contaminate the estimates since we set no constraint in the adjustment. The other potential reason could be that some other forces are not considered in the adjustment, for instance, the radiator effect in the Z direction. Also, there could be errors in satellite dimensions and total mass, especially for Block IIF satellites. Consequently, although the estimated optical parameters lack physical interpretation, the use of the estimated values can sufficiently model the SRP acceleration.

Table 2 The adjusted metadata for GPS Block IIA satellites (mass = 975 kg, sp denotes the solar panel, shape indicates the weighting between a cylindrical surface (1) and a flat surface (0))
Full size table
Table 3 Adjusted metadata for GPS Block IIR satellites (mass = 1100 kg, sp denotes the solar panel, the shape indicates the weighting between a cylindrical surface (1) and a flat surface (0))
Full size table
Table 4 Adjusted metadata for GPS Block IIF satellites (mass = 1555 kg, sp denotes the solar panel and shape indicates the weighting between a cylindrical surface (1) and a flat surface (0))
Full size table

In addition to the parameters in the tables, we have estimated solar panel rotation lag and solar sensor bias as well. The rotation lag of solar panels is about 1.5 deg for Block IIA satellites and about 0.8 deg for both Block IIR and IIF satellites. The sensor bias estimates differ clearly from satellite to satellite. All the Block IIA and IIR satellites have positive bias estimates, with a mean value of about 1.8 deg for IIA satellites and 0.7 deg for IIR satellites. However, Block IIF satellites have both positive and negative estimates, with a mean value of about 0.05 deg. We know that Block IIA satellites have a constant bias of + 0.5 deg implemented and we observe this value in our estimation. But we do not observe a clear − 0.5 deg bias that is mentioned in Sect. 2 for Block IIF satellites. The reason is not yet fully clear. If we look into the same Block IIF satellite (SVN 62, PRN 25) that has been analyzed by Dilssner (2010), we find a sensor bias estimate of − 0.3 deg. Since the solar panel rotation lag and the yaw bias effects can be well absorbed by the ECOM parameters (B sine and Y0), it is not fully necessary to consider them in the box-wing model.

GSPM model

The JPL GSPM model is a truncated Fourier expansion in elongation (Bar-Sever and Kuang 2004), as follows:

$$ {\mathbf{acc}}_{x} = s10^{ - 5} {{\left( {{{AU} \mathord{\left/ {\vphantom {{AU} r}} \right. \kern-\nulldelimiterspace} r}} \right)^{2} } \mathord{\left/ {\vphantom {{\left( {{{AU} \mathord{\left/ {\vphantom {{AU} r}} \right. \kern-\nulldelimiterspace} r}} \right)^{2} } m}} \right. \kern-\nulldelimiterspace} m}\left( {S_{ - } X_{1} \sin \varepsilon + S_{ - } X_{2} \sin 2\varepsilon + S_{ - } X_{3} \sin 3\varepsilon + S_{ - } X_{5} \sin 5\varepsilon + S_{ - } X_{7} \sin 7\varepsilon } \right) $$
(3)
$$ {\mathbf{acc}}_{y} = C_{ - } Y_{0} + 10^{ - 5} {{\left( {{{AU} \mathord{\left/ {\vphantom {{AU} r}} \right. \kern-\nulldelimiterspace} r}} \right)^{2} } \mathord{\left/ {\vphantom {{\left( {{{AU} \mathord{\left/ {\vphantom {{AU} r}} \right. \kern-\nulldelimiterspace} r}} \right)^{2} } m}} \right. \kern-\nulldelimiterspace} m}\left( {C_{ - } Y_{1} \cos \varepsilon + C_{ - } Y_{2} \cos 2\varepsilon } \right) $$
(4)
$$ {\mathbf{acc}}_{z} = s10^{ - 5} {{\left( {{{AU} \mathord{\left/ {\vphantom {{AU} r}} \right. \kern-\nulldelimiterspace} r}} \right)^{2} } \mathord{\left/ {\vphantom {{\left( {{{AU} \mathord{\left/ {\vphantom {{AU} r}} \right. \kern-\nulldelimiterspace} r}} \right)^{2} } m}} \right. \kern-\nulldelimiterspace} m}\left( {C_{ - } Z_{1} \cos \varepsilon + C_{ - } Z_{3} \cos 3\varepsilon + C_{ - } Z_{5} \cos 3\varepsilon } \right) $$
(5)

where \({\mathbf{acc}}_{x}\), \({\mathbf{acc}}_{y}\) and \({\mathbf{acc}}_{z}\) are accelerations in the satellite body frame XYZ (IGS definition), s is a dimensionless scaling factor (normally unity), AU the astronomical unit, r the distance between spacecraft and the Sun, m the spacecraft mass, and \(\varepsilon\) the elongation angle. Furthermore, \(C_{ - } Y_{0}\) is a constant Y bias and all the other terms are the empirical GSPM parameters.

The GSPM parameters are determined by JPL following three steps. First, select daily true orbit data from precise JPL orbit products and form 10-day orbit arcs. Satellites with large orbit fitting errors are excluded. Second, fit each 10-day-arc orbit data with the GSPM model and estimate the empirical parameters. Third, combine all satellite arc solutions into one set of GSPM parameters for Block IIA, IIR and IIF satellites, respectively. In the case of \(C_{ - } Y_{1}\) and \(S_{ - } X_{2}\), it is found that these two terms depend strongly on the \(\beta\) angle. A good fit between the observed variation and a function of \(\beta\) angle is suggested, as shown as follows:

$$ {\mathbf{F}}\left( \beta \right) = A + B\sin \beta + {C \mathord{\left/ {\vphantom {C {\sin \beta }}} \right. \kern-\nulldelimiterspace} {\sin \beta }} + D\cos \beta $$
(6)

where A, B, C, D are separate empirical parameters describing \(C_{ - } Y_{1}\) and \(S_{ - } X_{2}\) terms as a function of \(\beta\).

Different from fitting precise satellite orbits, we estimate daily GSPM parameters based on true GPS measurements of 6 years, similarly as for the box-wing parameters. Data inside eclipse seasons are not included since Block IIA satellites exhibit a highly nonlinear attitude behavior during eclipse seasons. We first generate daily solutions of 6 years. Then, we stack all the normal equations of 6 years, pre-eliminating all but the GSPM parameters. Finally, we obtain satellite-specific GSPM parameters and we combine all the solutions into Block-specific parameters. All the determined GSPM parameters are shown in Table 5.

Table 5 Estimates of the GSPM parameters for all GPS Blocks, in \(10^{ - 5}\) Newtons (note that all the parameters are given in the satellite body frame of the IGS definition)
Full size table

Validation

We name the box-wing model using optical properties from Rodriguez-Solano (2014) as the original box-wing (obw) model, while we call our adjusted box-wing model considering various physical effects on the enhanced box-wing (ebw) model. Qualities of these two box-wing models as well as the GSPM model are assessed by an orbit prediction over 5 days. The prediction is based on 1 day, 3 days, 5 days and 7 days of CODE final precise GPS orbit products. Only the constant ECOM Y0 parameter is employed in the orbit fitting. The predicted orbits of the 5th day are compared to the CODE final orbits of the same day. We compute RMS values of the 3D position differences of the orbit fitting and the orbit prediction. The prediction procedure is shifted day-by-day over a time period of 6 months of the year 2014 (day 001 to day 180).

Figures 8 and 9 show mean RMS values in orbit fitting and orbit prediction using different SRP models and different arc lengths. It is clear that the enhanced box-wing model performs much better than the original box-wing model, especially for Block IIA satellites. To interpret this, we would like to mention that we have observed a large correction of the estimated optical properties in the + X surface for Block IIA satellites. It seems that Rodriguez-Solano (2014) considers the area of the plume shield, a large flanged disk to protect the solar panels from apogee engine effluent, but does not consider its optical properties properly. Also, we have considered the radiator effect and the thermal radiation of the solar panel, which is another reason for the improvement. The GSPM parameters are determined using non-eclipsing data, and we try to validate them in the eclipse seasons as well. Results demonstrate that GSPM model shows similar performance as the enhanced box-wing model outside eclipsing seasons, whereas, in the eclipse season, the enhanced box-wing model performs better than the GSPM model for Block IIA and IIF satellites. As a consequence, the GSPM parameters determined in this section are only suitable for the non-eclipsing GPS satellites or for the eclipsing satellite with nominal satellite attitude, for instance, the Block IIR satellites. In the following section, we use the GSPM model only in non-eclipsing seasons.

Fig. 8
figure 8

RMS values of orbit fitting using different SRP models and different arc lengths. The top figure gives the result in the non-eclipse season, while the bottom figure gives the result in the eclipse season. To see the details, RMS values larger than 6 m are taken as 6 m

Full size image
Fig. 9
figure 9

RMS values of 5-day orbit prediction using different SRP models and different orbit fitting arc lengths. The top figure shows the result in the non-eclipse season, while the bottom figure shows the result in the eclipse season

Full size image

GPS satellite POD using individual SRP models

We take the original box-wing model, the enhanced box-wing model and the GSPM model as the a priori model, respectively, and jointly use them with ECOM and ECOM2 models in GPS satellite orbit determination. All the cases are listed in Table 6 where 7ECOM2 denotes the 7-parameter ECOM2 model without the fourth-order Fourier terms. The JPL GSPM model is, in fact, mostly used as an a priori model in a filter. Constant Y bias and SRP coefficient, as well as constrained time-varying SRP coefficients in satellite body-fixed X, Y and Z directions, are used to describe the SRP effect (http://acc.igs.org). We try to apply the GSPM model with ECOM models in the least-squares adjustment in this contribution.

Table 6 Different sets of SRP models. The GSPM model is only used in non-eclipse seasons
Full size table

We compute GPS satellite orbits using the same tracking stations and time periods as that in Sect. 4. The geodetic datum is defined by constraining station coordinates tightly to values in the IGS14 frame of the same station and epoch. Earth rotation parameters are estimated together with satellite orbit and clock parameters. Undifferenced phase ambiguities are fixed to integer values. Pseudo-stochastic pulses are not considered for the sake of assessing SRP models in a more dynamic mode. Orbit arc length is 3 days, and we extract the middle day of the 3-day-arc solution as the final daily solution. The ECOM parameters, orbit misclosures between two consecutive arcs, as well as orbit predictions of 24 h, are analyzed to assess the quality of the determined GPS satellite orbits.

Figures 10, 11 and 12 show the 5-parameter ECOM D0, Y0 and B0 estimates of three GPS satellites as a function of \(\beta\) angle. The results for different SRP models and different satellites are shifted by constant values for a clear display. It may be noted that the enhanced box-wing model and the GSPM model show more constant ECOM estimates than the other two cases except for the Block IIA satellite inside eclipse seasons. The reason is that the satellite attitude of IIA satellites inside shadows is not as accurate as that for IIR and IIF satellites. Radiator effects are thus not well compensated. RMS values of all the ECOM and ECOM2 parameters over 6 years for all the GPS satellites are shown in Table 7. The original box-wing model fails to model SRP accelerations in the ECOM B direction since we observe larger RMS values of obw + ECOM/7ECOM2/ECOM2 results than the others. The enhanced box-wing model and the GSPM model show similar RMS values; the mean improvement in all the parameters is about 30% compared to results without the a priori model.

Fig. 10
figure 10

ECOM D0 estimates of one Block IIA (SVN 10), one Block IIR (SVN 50) and one Block IIF (SVN 63) satellite. Mean values of each satellite are removed, and the results of different SRP models and different Blocks are shifted by constant values for better display

Full size image
Fig. 11
figure 11

ECOM Y0 estimates of one Block IIA (SVN 10), one Block IIR (SVN 50) and one Block IIF (SVN 63) satellite. Results of different SRP models and different Blocks are shifted by constant values for better display

Full size image
Fig. 12
figure 12

ECOM B0 estimates of one Block IIA (SVN 10), one Block IIR (SVN 50) and one Block IIF (SVN 63) satellite. Results of different SRP models and different Blocks are shifted by constant values for better display

Full size image
Table 7 RMS of the ECOM and ECOM2 parameters (unit: nm/s2); RMS values of the D0 parameter are the RMS around mean values
Full size table

Orbit misclosures at the day boundaries between consecutive arcs are traditionally used in orbit determination to assess the internal consistency of a specific processing scheme (Lutz et al. 2016). Figure 13 shows orbit misclosures of all the computed solutions as a function of \(\beta\) angle. Orbit misclosures of the ECOM model clearly show the largest values inside eclipsing seasons. The higher-order Fourier terms in the ECOM2 model can absorb part of the unmodeled forces and thus can reduce orbit errors inside the eclipse. The use of the original box-wing model has no positive effect on GPS satellite orbits. The GSPM model performs similarly as the enhanced box-wing model outside the eclipse. The enhanced box-wing model improves GPS satellite orbits in all the cases except for Block IIA satellites inside eclipse seasons. The reason is the same as we have discussed in Sect. 4 that the attitude inside eclipse seasons for Block IIA satellites is not as accurate as that for IIR and IIF satellites. RMS values of orbit misclosures for each solution are given in Table 8. The improvement in the enhanced box-wing model is clear. In particular, orbit misclosures of the eclipsing Block IIR and IIF satellites, as well as the non-eclipsing IIA satellites, is reduced by about 25%, 10% and 10%, respectively, for the ECOM model.

Fig. 13
figure 13

Orbit misclosures of all the GPS satellites over 6 years as a function of \(\beta\) using different SRP models (none denotes no a priori model, obw the original box-wing, ebw the enhanced box-wing and gspm the GSPM model

Full size image
Table 8 RMS of orbit misclosures for GPS satellite orbits using different SRP models (cm)
Full size table

Finally, we predict GPS satellite orbits over 24 h based on 3-day-arc orbit solutions (Duan et al. 2019a). The predicted orbits of different SRP models are compared to the CODE final orbits of the same time period. Table 9 shows RMS values of 3D orbit differences for each Block. Similar to orbit misclosures, the enhanced box-wing model shows improvement in all the cases except for that of Block IIA satellites inside eclipse seasons. GPS Block IIR satellites benefit clearly more from the enhanced box-wing model than Block IIA and IIF satellites. RMS values of 24-h predictions are reduced by about 10% and 20% for IIR satellites outside and inside eclipse seasons. The GSPM model exhibits no improvement for Block IIF satellites in the orbit prediction. The reason is that the GSPM parameters perform differently for each Block IIF satellite (Sakumura et al. 2017), and we simply adopt the mean values. Therefore, both orbit misclosures and 24-h orbit predictions demonstrate that the enhanced box-wing model is helpful for GPS satellite orbit determination for both ECOM and ECOM2 models.

Table 9 RMS of 3D orbit differences between 24-h predicted solutions and CODE final orbits for all the GPS satellites using different SRP models (cm)
Full size table

Summary and outlook

In this contribution, we propose an enhanced box-wing model considering various physical effects. We find that a yaw bias of + 0.5 deg for Block IIA satellites has an impact of about − 0.15 nm/s2 in the Y direction. By computing the ECOM Y0 estimates, we conclude that the radiator emission exists in the Y surface of Block IIA and IIR satellites. So, in addition to the traditional box-wing parameters, we adjust constant radiator effects in satellite body-fixed − X and Y directions. Other physical effects, such as thermal radiation and solar panel rotation lag, are considered as well. With all the adjusted box-wing parameters, we form an enhanced box-wing model.

To compare with the GSPM model, we determine Block-specific GSPM parameters using GPS measurements covering 6 years. Data inside eclipse seasons and outgassing periods are excluded. Performances of the original box-wing model, the enhanced box-wing model and the GSPM model are assessed by orbit prediction over 5 days for half a year. Both the enhanced box-wing model and the GSPM model show clear improvement compared to the original box-wing model. The GSPM model performs similarly as the enhanced box-wing model outside eclipse seasons, while it performs worse for Block IIA and IIF satellites inside eclipse seasons. As the GSPM parameters are determined from non-eclipsing data, they are only reliable for GPS satellites outside eclipse seasons. Since Block IIR satellites can maintain yaw nominal attitude inside shadows, the GSPM parameters show similar performance both inside and outside eclipse seasons.

The enhanced box-wing model and the GSPM model are then taken as a priori model and are jointly used with ECOM and ECOM2 models in GPS satellite orbit determination. The enhanced box-wing model performs similarly to the GSPM model outside eclipse seasons. Orbit misclosures of Block IIA satellites are improved by about 10% for the ECOM model. The enhanced box-wing model is quite efficient for the eclipsing satellites as well. Both orbit misclosures and 24-h orbit predictions demonstrate that the use of the enhanced box-wing model reduces RMS values by about 20% and 10% for Block IIR and IIF satellites inside eclipse seasons for the ECOM model.

In summary, the enhanced box-wing model is helpful for both ECOM and ECOM2 models. The ECOM2-related orbit misclosures and orbit predictions show clearly better performance than the 7ECOM2 and the ECOM model for results inside eclipsing seasons since higher-order Fourier terms can partly absorb some unmodeled forces. However, more analysis of station coordinates and geodetic parameters need to be done before concluding which SRP model is better.