1 Introduction

The GLONASS is operated by the Ministry of Defense of the Russian Federation. Since the launch of the first GLONASS satellite in 1982, a total of 24 satellites are now in operation (https://www.glonass-iac.ru). The first generation of GLONASS satellite was decommissioned in 2009, and the current GLONASS constellation consists of two types of spacecraft: GLONASS-M and GLONASS-K. Within the International GNSS Service (IGS), precise orbit and clock products of GLONASS satellites are routinely generated together with the US Global Positioning System (GPS) to support scientific and engineering applications (Dow et al. 2009; Johnston et al. 2017). Rather than GPS-only solutions, the combined GPS/GLONASS solutions provide more precise and reliable positioning results (Dach et al. 2009). Moreover, the initialization time in real-time kinematic applications is reduced significantly (Li et al. 2015). However, when assessing GLONASS-only solutions the Earth rotation parameters (ERPs) and geocenter coordinates demonstrate pronounced deviations, which were not present in the GPS-only solutions. It was shown by Meindl et al. (2013) and Rodriguez-Solano et al. (2014) that this is due to the imperfect SRP model of GLONASS satellites.

SRP is the dominant non-gravitational perturbation for Global Navigation Satellite System (GNSS) satellites. Milani et al. (1987) formulated the physical interaction between SRP and satellite surface by making use of attitude, total mass, dimensions, and optical properties of the satellites. However, since the metadata of GLONASS satellites are not yet officially published and we know very little about radiator and thermal radiation on the spacecraft, it is almost impossible to model SRP accelerations by using only analytical methods. In the absence of precise surface models, the ECOM model was developed in the early 1990s by Beutler et al. (1994), initially for the purpose of serving GPS satellites in yaw-steering mode. This empirical model was then modified and generalized by Springer et al. (1999) to serve as an a priori model for individual GPS satellites. In 2006, the CODE center calculated new coefficient sets for the model from Springer et al. (1999) by using CODE Final orbits from 2000 to 2006 and applied the same empirical model for GLONASS orbits as well (Dach et al. 2009). However, orbit predictions proved that the new GPS coefficients exhibited almost the same performance as before, while for GLONASS it was obviously best not to use the a priori model.

Different from GPS satellites, GLONASS-M satellites are cylindrical and have an elongated shape, and SRP acceleration varies periodically due to the asymmetric body dimensions, which cannot be sufficiently compensated by the classical ECOM model. Furthermore, GLONASS-M satellites assume a non-nominal attitude mode when entering into eclipse or experiencing noon turn maneuvers (Dilssner et al. 2011), and consequently, the effect of a radiator on a specific satellite surface might not be completely absorbed by the ECOM parameters. As an extension of the ECOM model, Arnold et al. (2015) developed the ECOM2 model, involving additional higher-order Fourier terms in the incident direction. They found that the spurious signals in the geocenter z-coordinate were significantly reduced for the GLONASS-only solutions compared to the ECOM results. Prange et al. (2017) applied ECOM and ECOM2 models for precise orbit determination of all the GNSS satellites. It was clear that the ECOM2 performed better than the ECOM model for Galileo and Quasi Zenith Satellite System (QZSS) satellite orbits, but the orbit misclosures of GLONASS orbits were minimally increased. All the GLONASS satellites inside the eclipse season were excluded from their analysis due to the particular attitude mode. Thus, although the higher-order Fourier coefficients can partly absorb unmodeled force errors, the ECOM2 model is not optimal for GLONASS orbit determination.

As proved by Steigenberger et al. (2015b), Montenbruck et al. (2015b, (2017), Zhao et al. (2018), Li et al. (2019), Duan et al. (2019b), an a priori box-wing model improves satellite orbits significantly if a satellite has an elongated shape or takes orbit normal mode. Based on this background, the goal of this paper is, therefore, to set up an improved a priori SRP model for GLONASS satellites. Without precise knowledge of the satellite metadata, optical parameters are adjusted as part of the precise orbit determination by making use of suitable dimensions and total mass of GLONASS satellites. The potential radiator and thermal radiation effects are taken into account. The adjusted optical parameters and radiator effects are introduced into an a priori box-wing model and jointly used with ECOM and ECOM2 models, respectively, in GLONASS POD.

2 Status of GLONASS POD results

We calculate GLONASS satellite orbits using 80 Multi-GNSS Experiment (MGEX) tracking stations; time span starts from day of year (doy) 200, 2017, to doy 60, 2018. Station-related parameters and ERPs are estimated together with satellite orbit and clock parameters. Pseudo-stochastic pulses are not considered for the sake of evaluating SRP models in a more dynamical mode, and we extract the middle day of the 3-day-arc solution as the final daily solution. GLONASS-M satellites employ a nominal yaw-steering mode outside the eclipse period, while take a specific attitude inside the eclipse and near the noon turn periods (Dilssner et al. 2011). GLONASS-K satellites are assumed to maintain yaw-steering mode all the time (Montenbruck et al. 2015a). Apart from the ECOM and ECOM2-only solutions, we set up an a priori box-wing model jointly used with the 5-parameter ECOM and 9-parameter ECOM2 model, respectively. Dimensions and total mass of satellites are taken from the IGSMAIL-5104 by Vladimir Mitrikas and Rodriguez-Solano (2014); optical parameters of satellite surfaces are adjusted following the same procedures as Duan et al. (2019b). Orbit processing settings are shown in Table 1.

Table 1 Settings of orbit determination
Full size table
Fig. 1
figure 1

The estimated ECOM parameters with (blue for GLONASS-M and yellow for GLONASS-K) and without (red for GLONASS-M and green for GLONASS-K) a priori model for GLONASS satellites. A mean value of \(150\,\hbox {nm}/\hbox {s}^2\) and \(104\,\hbox {nm}/\hbox {s}^2\) is added to the GLONASS-M and GLONASS-K \(D _{0}\) estimates without a priori model, respectively. The solid vertical lines indicate boundaries of eclipsing seasons

Full size image

The estimated ECOM parameters of GLONASS satellites as a function of \(\beta \) angle (the Sun elevation above the orbit plane) with and without a priori model are shown in Figure 1. It is clear that all the ECOM parameters perform differently inside the eclipse season. For the \(D _{0}\), \(B _{0}\) and Bc parameters, the use of the a priori box-wing model reduces the amplitude and the dependency on the \(\beta \) angle. However, the \(Y _{0}\) parameters exhibit symmetrical variations during the eclipse season, and the use of the a priori model has almost no effect. The potential reason can be that additional forces other than SRP contribute inside the shadow that cannot be compensated by the a priori box-wing or the ECOM model. When applying the ECOM2 model, the symmetrical variations of the \(Y _{0}\) parameters are reduced since the higher-order Fourier series can partly absorb unmodeled forces during eclipsing seasons.

We compute orbit differences at the day boundaries between consecutive arcs and define the orbit misclosure of two sets of orbits as

$$\begin{aligned} {m}_{i,i+1}=\sqrt{\frac{\sum _{s=1}^{n}|{{\textit{\textbf{r}}}}_{s,i+1}-{{\textit{\textbf{r}}}}_{s,i}|^2}{n}} \end{aligned}$$
(1)

where i and \(i+1\) refer to days, n the number of satellites available in both days, and \({{\textit{\textbf{r}}}}_{s,i+1}\) the orbit position of satellite s of day \(i+1\) (Lutz et al. 2016). Table 2 shows the mean values of orbit misclosures for individual SRP models inside and outside eclipse season. The ECOM model shows better performance than the ECOM2 model outside the eclipse season but is about two times worse inside the eclipse season. The a priori box-wing model results in almost no improvement for both ECOM models. Thus, forces inside the eclipse cannot be sufficiently accounted for by the ECOM model even if combined with a conventional box-wing model, but the ECOM2 model can absorb such forces fairly well.

Table 2 The mean values of orbit misclosures (unit : cm)
Full size table

All the GLONASS satellites are equipped with Laser Retroreflector Arrays (LRA) and are tracked by the International Laser Ranging Service (ILRS) network (Pearlman et al. 2002). The SLR range residuals are computed as differences between the SLR observations and the distances derived from the orbits. Table 3 summarizes the mean offsets and STD of the resulting one-way SLR residuals for all assessed solutions involving both eclipse and non-eclipse seasons. The STD value of the ECOM model is slightly better than of the ECOM2 model. Combined with an a priori box-wing model, the ECOM model results in the best SLR residuals.

Based on the status of GLONASS POD results, we find that the pure ECOM2 model is not optimal for GLONASS satellites, especially during the non-eclipse season. The conventional a priori box-wing model improves GLONASS orbits in the radial component, but fails to compensate for some forces inside the eclipse, which causes orbit errors in along-track and cross-track components.

Table 3 The mean offsets and STD of SLR residuals (unit : cm)
Full size table

3 New a priori SRP model

The GLONASS-M satellite body consists of two parts: a cylindrical container and boxlike antenna platform. As published by the GLONASS information center (www.glonass-iac.ru), GLONASS-M satellites have a width of 2.71 m and a length of 3.05 m. Rodriguez-Solano (2014) set up “shape” factors for the area of GLONASS-M satellites, which in fact represents the proportion of the cylindrical container and the boxlike antenna part. By making use of all the dimension information, we calculate a new area of the \(\pm z\) surface for GLONASS-M satellites. Note that the orientation of the satellite body frame in our study refers to the IGS definition (Montenbruck et al. 2015a).

3.1 Classical optical parameter adjustment

The physical interaction between the SRP and satellite solar panel is formulated as:

$$\begin{aligned} {{\textit{\textbf{f}}}}=-\frac{A}{M}\frac{S_{0}}{c}\cos \theta \bigg [(\alpha +\delta ){{\textit{\textbf{e}}}}_{D}+2\left( \frac{\delta }{3}+\rho \cos \theta \right) {{\textit{\textbf{e}}}}_{N})\bigg ] \end{aligned}$$
(2)

where A denotes the surface area, M the total mass of satellite, \(S_{0}\) the solar flux, c the vacuum velocity of light, and \(\alpha \), \(\delta \) and \(\rho \) the fractions of absorbed, diffusely scattered and specularly reflected photons. Furthermore, \({{\textit{\textbf{e}}}}_{D}\) denotes the Sun direction, \({{\textit{\textbf{e}}}}_{N}\) the surface normal vector, and \(\theta \) the angle between both vectors.

For the satellite body, surfaces are usually covered by multilayer insulation (MLI) blankets and the energy absorbed by the satellite surfaces is instantaneously reradiated back into space according to Lambert’s law.

$$\begin{aligned} {{\textit{\textbf{f}}}}=-\frac{A}{M}\frac{S_{0}}{c}\cos \theta \bigg [(\alpha +\delta )\left( {{\textit{\textbf{e}}}}_{D}+\frac{2}{3}{{\textit{\textbf{e}}}}_{N}\right) +2\rho \cos \theta {{\textit{\textbf{e}}}}_{N}\bigg ]. \end{aligned}$$
(3)

In case of a cylindrical satellite surface, the acceleration is formulated by Fliegel et al. (1992) as:

$$\begin{aligned} {{\textit{\textbf{f}}}}=-\frac{A}{M}\frac{S_{0}}{c}\cos \theta \bigg [(\alpha +\delta )\left( {{\textit{\textbf{e}}}}_{D}+\frac{\pi }{6}{{\textit{\textbf{e}}}}_{N}\right) +\frac{4}{3}\rho \cos \theta {{\textit{\textbf{e}}}}_{N}\bigg ]. \end{aligned}$$
(4)

The GLONASS-M satellite body contains both cylindrical and flat surfaces, and consequently, Rodriguez-Solano et al. (2012a) combined Eqs. (3) and (4) considering shape factor “s”.

$$\begin{aligned} {{\textit{\textbf{f}}}}=&-\frac{A}{M}\frac{S_{0}}{c}\cos \theta \left[ (\alpha +\delta )\left( {{\textit{\textbf{e}}}}_{D}+\left( \frac{\pi }{6}s+\frac{2}{3}(1-s)\right) {{\textit{\textbf{e}}}}_{N}\right) \right. \nonumber \\&\left. +\left( \frac{4}{3}s+2(1-s)\right) \rho \cos \theta {{\textit{\textbf{e}}}}_{N}\right] . \end{aligned}$$
(5)

In our estimation, Eq. (5) is applied in a box-wing model to study the satellite body surfaces, while Eq. (2) is applied for the solar panel. Since GLONASS satellites take yaw-steering mode in most of their mission, only optical parameters in \(+x\) and \(\pm z\) panels are adjusted. For each satellite surface, there are two adjusted optical parameters: \(\alpha +\delta \) and \(\rho \). In order to compensate for the direct radiation and the misalignment of the solar panel, we estimate additionally a scaling factor, y-bias and a rotation lag of the solar panel (Rodriguez-Solano et al. 2012a, b). We use the same data as in Sect. 2; days with maneuvers are excluded from the adjustment. Station coordinates, troposphere delays and receiver clock offsets are fixed to GPS precise point positioning (PPP) solutions (Steigenberger et al. 2013, 2015a; Selmke and Hugentobler 2017). Dimensions and total mass of the satellites are fixed to the given values; the initial optical parameters from Rodriguez-Solano (2014) are taken as a priori values in the least squares adjustment. Daily normal equations are then stacked after pre-eliminating all but the optical parameters. Finally, we take the averaged values of all the GLONASS-M estimates and there is no need to do a second iteration. There is only one GLONASS-K satellite (R09) continuously available during the experimental period, and thus the optical parameters of GLONASS-K satellite are the estimates of one satellite.

Fig. 2
figure 2

Daily acceleration estimates in \(-x\) direction and formal uncertainties for one GLONASS-M (using dedicated attitude mode presented by Dilssner et al. (2011)) and one GLONASS-K satellite (using additionally nominal yaw-steering mode) inside the shadow

Full size image
Fig. 3
figure 3

Daily radiator parameter estimate in \(-x\) direction for all GLONASS satellites; the solid vertical lines indicate boundaries of eclipsing seasons

Full size image
Fig. 4
figure 4

Simulated accelerations of a GLONASS-M satellite (\(\beta =5\) deg) by the \(+x\) surface in the x direction with and without self-shadowing

Full size image
Table 4 Different sets of GLONASS metadata
Full size table
Fig. 5
figure 5

The mean RMS values of orbit differences between the predictions of the 7th-day and the precise CODE orbits for each GLONASS satellite by using individual types of metadata; a, b, c indicate individual orbit planes with averaged \(\beta \) angle of \(43^{\circ }\), \(23^{\circ }\) and \(-48^{\circ }\), respectively; GLONASS-K is R09

Full size image

3.2 Thermal radiation and radiator effects

In the previous adjustment of optical parameters, for the sake of analytical simplicity, temperatures of both sides of the solar panels are assumed to be the same. However, in the real case, temperature of the front side should be higher than the back side, but in general we know little about the details of the temperature difference. As the density of Earth’s atmosphere in space at GNSS’s altitude is extremely low, there is almost no heat transfer by heat convection, the only transfer mechanism between the solar panel and the environment is thermal radiation (Walter 2012). Since the incident vector is always perpendicular to the solar panel in the yaw-steering mode, the total effect of thermal radiation can be assumed to be in the same direction as the effect of the photons reflected at the solar panel. Therefore, in the optical parameter adjustment, we fix \(\alpha \) and \(\delta \) of the solar panel on the a priori values and estimate a correction of the parameter \(\rho \) of the solar panel.

Spacecraft needs radiators to get rid of excess heat, and in most cases, their location is selected to minimize incident radiation. Thus, potential radiators of GLONASS satellites can be expected on the \(-x\) and \(\pm y\) surfaces since these are not illuminated in the yaw-steering mode. Results of the \(Y _{0}\) parameters in Sect. 2 demonstrate that y-bias of all the GLONASS satellites is constant outside the eclipse seasons with an acceleration of about \(0.2\,\hbox {nm}/\hbox {s}^2\). Therefore, radiator effects in \(\pm y\) surfaces are either balanced or are very small. However, \(Y _{0}\) parameters show symmetrical variations during the eclipse season with an amplitude of about \(1\,\hbox {nm}/\hbox {s}^2\) for GLONASS-M satellites and \(0.5\,\hbox {nm}/\hbox {s}^2\) for the single GLONASS-K satellite. The reason might be that radiators on the \(-x\) surface contribute as well when satellites are in the shadow. Since the attitude of GLONASS satellites inside the shadow is far from the yaw-steering mode, the effect of radiators on the \(-x\) surface causes an acceleration in the nominal y-direction, which is wrongly absorbed by the ECOM \(Y _{0}\) parameter. Nevertheless, the higher-order Fourier terms in \({{\textit{\textbf{e}}}}_{D}\) can partly absorb such radiator effects. This is why the ECOM2 model performs better than the ECOM model inside the eclipse seasons. Therefore, in our optical parameter adjustment we add one additional constant radiator parameter in the \(-x\) direction. However, the attitude control law of GLONASS-K satellites inside the shadow is unknown since we cannot apply the same reverse-kinematic-positioning-like approach as for GLONASS-M satellites due to the fact that the horizontal satellite antenna offset of GLONASS-K satellites is zero.

Table 5 The adjusted metadata of GLONASS-M satellites (mass = 1415 kg, sp denotes solar panel, shape indicates the weighting between a cylindrical surface (1) and a flat surface (0))
Full size table
Table 6 The adjusted metadata of GLONASS-K satellites (mass = 935 kg, sp denotes solar panel)
Full size table
Fig. 6
figure 6

The estimated ECOM parameters with (blue for GLONASS-M and yellow for GLONASS-K) and without (red for GLONASS-M and green for GLONASS-K) new a priori model for GLONASS satellites. A mean value of \(150\,\hbox {nm}/\hbox {s}^2\) and \(104\,\hbox {nm}/\hbox {s}^2\) is added to the GLONASS-M and GLONASS-K \(D _{0}\) estimates without a priori model, respectively. The solid vertical lines indicate boundaries of eclipsing seasons

Full size image

In this contribution, we present a new approach allowing to confirm that GLONASS-K satellites take a similar attitude mode as GLONASS-M satellites inside the shadow. For this, we set up a constant acceleration parameter in the \(-x\) direction inside the shadow and jointly estimate it together with the ECOM parameters as part of POD using both nominal and dedicated GLONASS-M attitudes. Since satellite observations base on the “true” satellite attitude, we can then evaluate the performance of the estimated accelerations inside the shadow. Figure 2 shows the daily estimates and uncertainties of one GLONASS-M satellite and one GLONASS-K satellite in the eclipse season. The GLONASS-M satellite takes the attitude mode presented by Dilssner et al. (2011), while GLONASS-K satellite assumes in addition the nominal attitude mode. At the edges of eclipse season, the shadow time period is short and the accelerations are not well determined due to the low number of observations. The dedicated GLONASS-M attitude is confirmed for GLONASS-M satellites; the variation of daily estimates is small. When applying the same attitude mode to the GLONASS-K satellite, we observe similar performance as for the GLONASS-M satellite. However, the result for the GLONASS-K satellite using nominal attitude shows clearly a lager scatter, and the formal uncertainty is also in general two to three times larger. Therefore, we can conclude that GLONASS-K satellites employ a similar attitude control law as GLONASS-M satellites inside the shadow. In the following optical parameter adjustment, the dedicated GLONASS-M satellite attitude mode inside the eclipse is applied for GLONASS-K satellite as well. The radiator effect in \(-x\) direction is considered as a daily parameter for each satellite both inside and outside the shadow and is jointly estimated together with other optical parameters. To investigate the daily variation of the radiator estimates, we stack optical parameters in the normal equations but provide daily estimates of radiator acceleration, as shown in Fig. 3. The estimates show a slightly “curved” variation as a function of the \(\beta \) angle. The formal errors are getting large if the \(\beta \) angle is high since the radiator parameter is correlated with the direct scaling factor parameter. In our adjustment, we assume that the radiator parameter of one satellite does not vary with time and we stack them in the normal equations together with all the optical parameters. Furthermore, the GLONASS-M satellites take the averaged radiator estimates of all the GLONASS-M satellites, while the value for GLONASS-K is based on a single satellite solution. As a result, the radiator and optical parameters are used to set up an improved a priori box-wing model.

In addition, the antenna platform of GLONASS-M satellites is larger than the diameter of the cylindrical container. Hence, if the \(\beta \) angle is smaller than a threshold there would be a shadow caused by the antenna onto the \(+x\) surface of the cylindrical container. By making use of all the dimension information, we assume that the diameter of the cylinder is 1.63 m, the additional part of the antenna panel above the cylindrical container is 0.5 m, and the length of the cylindrical container is 2.03 m. Then, we can calculate the effects of self-shadowing for absorption, emission and reflection separately. Figure 4 shows the accelerations caused by the \(+x\) surface in the x direction with and without considering shadow effect for one GLONASS-M satellite with a \(\beta \) angle of \(5^{\circ }\). It is observed that the effect of shadow can be as large as \(1\,\hbox {nm}/\hbox {s}^2\) in the x direction.

Fig. 7
figure 7

The estimated ECOM2 parameters with (blue for GLONASS-M and yellow for GLONASS-K) and without (red for GLONASS-M and green for GLONASS-K) new a priori model for GLONASS satellites. A mean value of \(150\,\hbox {nm}/\hbox {s}^2\) and \(104\,\hbox {nm}/\hbox {s}^2\) is added to the GLONASS-M and GLONASS-K \(D _{0}\) estimates without a priori model, respectively. The solid vertical lines indicate boundaries of eclipsing seasons

Full size image

3.3 Newly adjusted optical parameters

Considering the above-discussed factors, we have four sets of metadata (dimensions, total mass and optical properties) for GLONASS satellites, as shown in Table 4. The first set is consistent with Rodriguez-Solano (2014), the second is the same as the first except for the new area of the \(\pm z\) surface, the third set is the same as the second but with adjusted optical parameters considering the thermal radiation and radiator effects, and the fourth set is the same as the third but considering additionally the self-shadowing effect on the \(+x\) surface. Dimensions of GLONASS-K satellite are equal in all cases, and there is no self-shadowing; thus, set 1 is the same as set 2 and set 3 is the same as set 4 for GLONASS-K satellites.

The quality of the adjusted parameters is assessed by an orbit prediction over one week. We fit CODE orbits of one day to generate initial conditions without considering any empirical parameters. Accelerations of SRP are calculated based on individual sets of metadata. The predicted orbits of the 7th day are compared to the precise CODE orbit products of the same day. This processing procedure is shifted day by day over a time interval of two weeks of the year 2017. GLONASS satellites in the maneuver or eclipse periods are excluded from the analysis. Figure 5 shows the quality of the predicted orbits in radial, along-track and cross-track components by using individual sets of metadata. It is clear that set 2 is better than set 1, which indicates that the new area of the \(\pm z\) surface is more reasonable than the currently known value. Set 3 and set 4 are in general 2 to 4 times better than set 2, but the difference between set 3 and set 4 is quite small. To further assess the self-shadowing effect, we do the same experiment for periods where one orbit plane has a \(\beta \) angle close to zero. We find that there is an improvement of about 10% in the along-track direction by considering the self-shadowing effect. In fact, we have also applied both set 3 and set 4 metadata into precise orbit determination, but the results are quite similar. Therefore, although the self-shadowing effect of \(+x\) surface is correctly modeled the impact on the precise orbit determination is not significant. Since the calculation of the self-shadowing effect in set 4 is complicated and no accurate official dimensions of the satellites are published, we recommend to use set 3 metadata for precise orbit determination of GLONASS satellites.

Table 7 RMS of the ECOM and ECOM2 parameters (unit : \(\hbox {nm}/\hbox {s}^2\) ). Note that the RMS values of the D0 parameter are the RMS around the mean value
Full size table
Table 8 The mean values of orbit misclosures by using the new box-wing model; 7-ECOM2 denotes 7-parameter ECOM2 and ECOM2 represents 9-parameter ECOM2 (unit : cm)
Full size table
Table 9 The mean RMS values of orbit differences between 24-h predicted and adjusted orbits; R means radial component, A denotes along-track component, C means cross-track component, 7-ECOM2 denotes 7-parameter ECOM2, and ECOM2 represents 9-parameter ECOM2 (unit : cm)
Full size table

The corresponding estimated metadata are given in Tables 5 and 6. As the \(-x\) panel is not illuminated, optical properties needed for albedo are assumed to be the same as for \(+x\). All the formal uncertainties are reasonably small. A radiator acceleration of about \(-1\,\hbox {nm}/\hbox {s}^2\) is observed for GLONASS-M satellite, which indicates roughly 500 W emission power. Unphysical negative optical properties are found for some surfaces. The first potential reason is that parameter \(\alpha +\delta \) and \(\rho \) of a satellite surface are correlated, and systematic errors might contaminate the estimates since we set no constraint in the adjustment. To reduce the correlation effect, we have attempted to constrain all the negative values to zero and got a new set of estimates. However, when assessing the 7-day predicted orbits we find that the constrained estimates lead to less good results than the estimates without any constraints. Therefore, apart from the correlation effect, there could be some other unmodeled forces affecting our adjustment, for instance radiators on the \(\pm z\) surface. Also, the sum of \(\alpha +\delta \) and \(\rho \) for each satellite surface is not equal to the physical condition 1. Nevertheless, the summed values of the \(\pm z\) surfaces are almost the same, which indicates an identical scaling factor. Consequently, the deviation of the summed values from 1 might be due to errors in the dimensions or the total mass of the satellites. Thus, although the estimated optical parameters lack physical interpretation, the use of the radiator and optical parameters seems to sufficiently model the SRP accelerations.

4 New GLONASS POD results

With the new a priori box-wing model, we estimate GLONASS satellite orbits using the same data and processing settings as in Sect. 2. The empirical parameters, orbit misclosures between two consecutive arcs, orbit predictions and SLR residuals are analyzed to assess the quality of GLONASS orbits.

4.1 The empirical parameters

If the a priori box-wing model is perfect, the ECOM estimates should be zero. Starting from this fact, we evaluate the effect of the introduced new box-wing model. The ECOM and ECOM2 parameters with and without the new a priori box-wing model are shown in Figs. 6 and 7, respectively. A direct acceleration of \(-150\,\hbox {nm}/\hbox {s}^2\) and \(-104\,\hbox {nm}/\hbox {s}^2\) is subtracted from GLONASS-M and GLONASS-K \(D _{0}\) parameters estimated without a priori model.

Table 10 SLR residuals by using the new box-wing model; 7-ECOM2 denotes 7-parameter ECOM2 and ECOM2 represents 9-parameter ECOM2 (unit : cm)
Full size table
Fig. 8
figure 8

SLR residuals of GLONASS satellite orbits as a function of \(\beta \) and \(\Delta u\) by using ECOM and BW+ECOM models

Full size image
Fig. 9
figure 9

SLR residuals of GLONASS satellite orbits as a function of \(\beta \) and \(\Delta u\) by using 7-ECOM2 and BW+7-ECOM2 models

Full size image

The ECOM and ECOM2 \(D _{0}\) parameters without a priori model show a clear dependency on the \(\beta \) angle. When considering the new a priori model, the dependency is significantly reduced, and \(D _{0}\) parameters are constant also during the eclipse period. Different from the symmetrical variations of ECOM \(Y _{0}\) parameter in Fig. 1, the use of the new box-wing model reduces such variations clearly in the eclipse season for both GLONASS-M and GLONASS-K satellites. The maximum acceleration drops from as large as about \(1\,\hbox {nm}/\hbox {s}^2\) to less than \(0.5\,\hbox {nm}/\hbox {s}^2\). The reason is the partially modeled effect of the radiator on the \(-x\) surface during the eclipse. The ECOM2 \(Y _{0}\) parameters show smaller variations than the ECOM \(Y _{0}\) during the eclipse season since the high-order Fourier terms can partly absorb errors in the force model. For all the other empirical parameters, the use of the new box-wing model reduces the amplitude and the dependency on the \(\beta \) angle. Table 7 shows the RMS values of ECOM and ECOM2 parameter series over the whole time period. In general, the ECOM model benefits more than the ECOM2 model, with an improvement of about a factor of two, which indicates that the adjusted optical parameters, radiator and thermal radiation effects are reasonable. Furthermore, we find that the uncertainties of ECOM2 D4C/D4S term are larger than the other ECOM2 estimates, which indicates that the ECOM2 solutions might be overparameterized. Also, as presented by Dach et al. (2016) the fourth-order terms are not recommended. In the following subsection, we show 7-parameter ECOM2 results excluding D4C and D4S terms as well.

Fig. 10
figure 10

SLR residuals of GLONASS satellite orbits as a function of \(\beta \) and \(\Delta u\) by using ECOM2 and BW+ECOM2 models

Full size image
Fig. 11
figure 11

SLR residuals of GLONASS satellite orbits as a function of elongation by using ECOM and BW+ECOM models

Full size image
Fig. 12
figure 12

SLR residuals of GLONASS satellite orbits as a function of elongation by using 7-ECOM2 and BW+7-ECOM2 models

Full size image
Fig. 13
figure 13

SLR residuals of GLONASS satellite orbits as a function of elongation by using ECOM2 and BW+ECOM2 models

Full size image

4.2 Orbit validation

Similar as in Sect. 2, we compute orbit misclosures at the day boundaries between consecutive arcs. Table 8 shows the mean values of orbit misclosure by using the new box-wing model inside and outside the eclipse season. Different from the conventional box-wing model, the newly employed a priori box-wing model reduces orbit misclosures in all cases for both GLONASS-M and GLONASS-K satellites. In particular, the improvement in the ECOM model is more than 30% in the eclipse seasons, and in general the BW+ECOM model results in the best orbit misclosures.

Then, we predict GLONASS satellite orbits over 24 hours based on 3-day-arc orbit solutions using different SRP models (Duan et al. 2019a). The predicted orbits are compared to the corresponding adjusted orbits for the same period. Table 9 shows the mean RMS values of orbit differences. The new a priori box-wing model reduces the RMS values for both the ECOM and ECOM2 models inside and outside the eclipse season. In particular, the improvement in the ECOM model for GLONASS-M satellites is about a factor of two in the eclipse seasons. In general, the BW+ECOM model shows the best 24-h predictions outside the eclipse seasons, while the BW+ECOM2 model results in better predictions in along- and cross-track directions inside the eclipse seasons. The GLONASS-K satellite shows less improvement than GLONASS-M satellites in both orbit misclosures and predictions. The reason is that the shape of GONASS-K satellite is less elongated and the radiator effect is smaller.

Finally, the adjusted GLONASS satellite orbits are validated by SLR observations. Figures 8, 9 and 10 illustrate the SLR residuals as a function of \(\beta \) and \(\Delta u\), where \(\Delta u\) is the difference between the satellite’s and the Sun’s argument of latitude. For the ECOM model, the maximum positive bias is found when \(\Delta u\) is close to \(180^{\circ }\), while the maximum negative bias is found when \(\Delta u\) is close to zero. Combining the ECOM model with the new a priori box-wing model greatly reduces the spurious pattern of the SLR residuals. For the ECOM2 model, such a pattern is not as obvious as for the ECOM model, but still we observe a maximum negative bias when \(\Delta u\) is close to zero. When combining with the new box-wing model, such a pattern disappears. Figures 11, 12 and 13 show the dependency of SLR residuals on the elongation (the angle between Earth and Sun as seen by the satellite). It is clear that the slopes for both ECOM and ECOM2 models reduce to almost zero. As a result, the use of the new box-wing model makes the estimated orbits almost unaffected by artifacts related to SRP modeling deficiencies.

Table 10 summarizes the mean offsets and STD values of the SLR residuals for all assessed solutions. The STD value of the ECOM model and 7-parameter ECOM2 model is slightly better than for the ECOM2 model. The new a priori box-wing model reduces the mean and STD values of SLR residuals for all the ECOM models, and GLONASS-M satellites gain more than the GLONASS-K satellite. In summary, the BW+ECOM model results in the best SLR residuals, and the mean offset and the STD value reduce by about 50% and 20%, respectively, with respect to the orbit products without the a priori model.

5 Summary and conclusions

GLONASS satellite orbit products are routinely generated within the IGS. In the current state of the art, either the ECOM or ECOM2 model is employed with stochastic pulses at the middle of a day, and GLONASS-M satellites are known to take a specific attitude mode during eclipse seasons, while GLONASS-K satellites are assumed to take yaw-steering mode all the times. In this contribution, we first perform the conventional processing and confirm that GLONASS-K satellites follow similar attitude mode as GLONASS-M satellites inside the eclipse. Different from the conventional box-wing model, potential radiators on the \(-x\) panel and thermal radiation of the solar panel are taken into account in the optical parameter adjustment as part of the precise orbit determination. From the adjusted parameters, a new box-wing model is conducted that is jointly used with the ECOM and ECOM2 model, respectively.

Considering the geometry of GLONASS-M satellites, we have developed a self-shadowing model, allowing us to calculate the effect of the shaded area on the \(+x\) cylindrical surface caused by the antenna part. We have proved that the improvement is about 10% in the along-track component for the eclipsing satellites when predicting satellite orbits of 7 days without any empirical parameters. However, the effect is not significant in the satellite orbit determination and thus can be ignored.

Comparing with the pure ECOM models and the conventional box-wing model, the advantage of the new box-wing model is clear. The amplitude and the dependency of all the empirical parameters on the \(\beta \) angle are greatly reduced, which indicates that the new box-wing model represents GLONASS satellites correctly. Orbit misclosures and RMS values of 24-h predicted orbits are all improved by using the new box-wing model. In particular, the improvement is more than 30% for the ECOM model during the eclipsing seasons. Finally, satellite orbits are assessed by SLR observations. The use of the new box-wing model reduces the spurious pattern of the SLR residuals as a function of \(\beta \) and \(\Delta u\) significantly, and the linear dependency of the SLR residuals on the elongation drops from as large as \(-0.760\) mm/deg to almost zero for both ECOM and ECOM2 models. The BW+ECOM model results in the best SLR residuals, with a mean and STD value of 1.1 and 4.4 cm for GLONASS-M satellites and 2.6 and 4.6 cm for the GLONASS-K satellite. In general, the GLONASS-K satellite benefits less than GLONASS-M satellites from the new a priori box-wing model due to the less elongated body shape and smaller radiator effect.