1 Introduction

1.1 The European navigation satellite system

The European navigation satellite system Galileo is being developed by the European Space Agency (ESA) under the partnership of the European Union and is the first, fully civilian navigation system. After a successful phase of the Galileo In-Orbit Validation Elements (GIOVE-A and -B, Montenbruck et al. 2006; Steigenberger et al. 2011), four In-Orbit Validation (IOV) spacecraft have been launched. To this day, three satellites are available in the operational constellation because one of the IOV satellites, GAL-104, experienced a sudden power loss which resulted in the permanent failure of E5 and E6 signal transmission; therefore, GAL-104 is available only for single-frequency navigation (Steigenberger and Montenbruck 2017).

Table 1 Characteristics of the Galileo satellites and their orbit features
Full size table

The operational phase of the Galileo system is being introduced since the launch of the first pair of the Full Operational Capability (FOC) satellites. However, the first two satellites have been accidentally placed on highly eccentric orbits (Sośnica et al. 2018). Despite the fact that the two satellites are not fully suitable for navigation, they may serve for the investigation of gravitational redshift (Delva et al. 2015; Herrmann et al. 2018), as well as for precise positioning (Paziewski et al. 2018). The following launches resulted in the increase in the number of active Galileo satellites allowing the European GNSS Agency (GSA) for announcing the operational capability in December 2016 when the part of operational services have been activated. Currently, the Galileo constellation consists of 26 spacecraft, i.e., four IOV satellites and 22 FOC satellites (out of which two fly at highly elliptical orbits and one, GAL-204, has been removed from the operational services due to the clock issues). All the Galileo satellites are equipped with Laser Retroreflector Arrays (LRA) for Satellite Laser Ranging (SLR) and are tracked by the network of the International Laser Ranging Service (ILRS, Pearlman et al. 2019). The current status of the Galileo constellation together with the Galileo orbit characteristics is shown in Table 1.

The Galileo satellites follow in principle the yaw-steering mode in which the bus of the satellite is fixed in a way it keeps the Z-axis toward the Earth center, thus illuminates the Earth with its navigation signal. The Y-axis is perpendicular to the Sun and the X-axis points toward the clock panel which is directed into deep space. In order to maintain the nominal attitude, it is necessary to turn (‘yaw’) about its Z-axis while rotating its solar panels around the Y-axis. In this study, we use the IGS nomenclature (Montenbruck et al. 2015b) where the Galileo body frame is inverted w.r.t the manufacturer-specific frame, i.e., as shown in Fig. 1. The non-gravitational forces are usually considered in the Sun-satellite-Earth (SSE) reference frame. The Sun-induced forces are dependent on the Sun elevation above the orbital plane, that is \(\beta \) angle, as well as the argument of latitude of the satellite w.r.t the argument of latitude of the Sun (\(\varDelta \)\(u\)), and the elongation angle \(\varepsilon \) (see Fig. 1).

Fig. 1
figure 1

Sun-satellite-Earth reference frame illustrating the angular elevation of the Sun above the orbital plane (\(\beta \)), the argument of latitude of the satellite with respect to the argument of latitude of the Sun (\(\varDelta u\)), and the elongation angle \(\epsilon \). The X, Y, and Z directions are consistent with IGS (Montenbruck et al. 2015b)

Full size image

1.2 Non-gravitational forces acting on GNSS satellites

According to Marquis and Krier (2000), non-gravitational perturbing forces acting on GNSS satellites result primarily from: direct solar radiation pressure (SRP), radiation of thermal blankets, thermal radiation from the satellite radiators, solar panels thermal radiation, thermal radiation of excess solar array power (shunt), and the Earth radiation pressure (ERP) which can be divided into the reflected Earth’s surface solar radiation (albedo), and Earth’s infrared radiation (IR). Moreover, GNSS satellites transmit continuously the radio navigation signals. The signal transmission causes systematic accelerations continuously acting on the satellite which can be identified as the navigation antenna thrust (Steigenberger et al. 2017).

The visible contribution of the Sun comprises more than 95% of the total force. The second highest is albedo reaching 2.5%. The solar array and shunt thermal radiation forces comprise less than 1% of the impact. The magnitudes of the two last forces are approximately the same; however, they act in opposite directions, thus their impact cancels each other (Marquis and Krier 2000). The aforementioned analyses were performed for the GPS satellites, however, according to the Galileo metadata, the radiators mounted on the +/\(-Y\) panels of the Galileo-FOC satellites are asymmetrical which may cause systematic effects in their orbit modeling.

Coping with SRP, together with the albedo, IR, satellite thermal effects and the impact of the navigation antenna thrust is crucial also in terms of the Galileo precise orbit determination. Direct SRP acting on the GNSS satellites can be handled threefold (Ziebart 2004): using empirical orbit models such as an Empirical CODE Orbit Model (ECOM2, Arnold et al. 2015), analytical models, e.g., using ray-tracing technique (Li et al. 2018) or hybrid models such as an a priori cuboid model proposed by Montenbruck et al. (2015a).

Empirical models typically decompose the accelerations acting on the GNSS satellites in three directions in the SSE frame (Fig. 1), i.e., D—pointing from the satellite toward the Sun, Y—along the solar panel rotation axis, and B—perpendicular to D and Y axes, completing the right-handed orthogonal frame. The original ECOM model was proposed by Beutler et al. (1994). However, with the emergence of new GNSS constellations, the ECOM model became obsolete due to new shapes of satellite buses. Therefore, a modified ECOM2 has been introduced by Arnold et al. (2015). The current version of the ECOM2 considers the constant accelerations in all directions (D0, Y0, B0), even sine and cosine terms in direction D, currently limited to D2C and D2S, and odd terms in the direction B, that is B1C and B1S. Empirical models can also absorb some of the accelerations which result from both albedo and IR to the extent which stems from the set of the estimated parameters. However, the accelerations resulting from ERP act on the GNSS satellites in different periods as compared to SRP. Therefore, empirical models are not capable of describing the physical interaction between the perturbing forces and the satellite surfaces.

The analytical approach can be applied with the adaptation of solely physical properties of the satellites using, for instance, the ray-tracing technique (Li et al. 2018) or in a hybrid way, with the additional estimation of the empirical parameters. Physical properties allow for the composition of the ‘box-wing’ model which simplifies spacecraft to the satellite bus (‘box’) and solar panels (‘wing’). Such models have been proposed for the GPS satellites by Rodriguez-Solano et al. (2012a) and for the old type of GLONASS satellites (Ziebart and Dare 2001). However, for both systems there are no official optical parameters available; hence, the box-wing model parameters have to be adjusted during the observation processing. When the optical properties are unknown, it is possible to formulate the approximate cuboid box-wing model such as the one provided by Montenbruck et al. (2015a)

The feasible usage of the analytical or semi-analytical models for the Galileo satellites is possible since late 2017 when GSA released the physical properties for the Galileo constellation.Footnote 1 Tests of the optical properties were performed by Duan et al. (2018), whereas the overview of the precise orbit determination strategies based on the Galileo metadata are studied by Bury et al. (2019) and Li et al. (2019). Based on the optical properties, it is possible to evaluate the accelerations resulting from SRP, albedo, and IR and take them into account for the Galileo precise orbit determination (Rodriguez-Solano et al. 2012b). Prange et al. (2017a) indicated that median SLR residuals decrease by 18 mm when ERP modeling is applied. However, none of the studies provides a complex analysis of the characteristics and the impact of the particular non-gravitational perturbing forces on the Galileo satellite orbits.

The latest study on the antenna thrust impact was performed by Steigenberger et al. (2017). Their analyses, based on the antenna power value measurement, indicated that for the Galileo satellites, the negligence of the antenna thrust introduces a significant shift in the radial direction at the level of 27.0 and 24.7 mm for the Galileo-FOC satellite on nominal and eccentric orbits, respectively. In terms of the Galileo-IOV satellites, the systematic shift is at the level of 13.9 and 9.8 mm for the operational IOV satellites, and for the E20 satellite whose power was reduced due to power supply issues. A similar analysis was performed by Prange et al. (2017a), who evaluated the impact of the antenna thrust on the Galileo orbit as of 1 cm per 100 W power.

The majority of the perturbing forces resulting from SRP and ERP can be absorbed owing to the application of the box-wing model. This brings us closer to the 1-cm level Galileo orbit which is crucial in the light of the International GNSS Service (IGS, Johnston et al. 2017) Multi-GNSS Pilot Project (MGEX, Montenbruck et al. 2017) whose main goal is the development of the orbit strategy for all new emerging navigation systems.

1.3 Goal of this study

The goal of this study is to evaluate the impact of selected, the most crucial for the precise orbit determination, non-gravitational forces acting on the Galileo satellites. The impact of the direct SRP, albedo, and Earth IR is assessed based on the box-wing model which is composed with the usage of the official metadata released by GSA. Additionally, we check the influence of the navigation antenna thrust on the Galileo satellite orbits. Employing different orbit determination strategies, we evaluate the box-wing model efficiency based on the analysis of the empirical parameters estimated with and without the box-wing model. Finally, we indicate the best performing strategy for the Galileo orbit determination applying the aforementioned force model.

The methodology of the orbit determination is described in Sect. 2. The magnitude and periodicity of accelerations resulting from SRP, albedo, and IR is described in Sect. 3. Moreover, Sect. 3 provides information about the sensitivity of the Galileo satellites on both Y- and B-biases. Precise orbit determination of the Galileo satellites is described in Sect. 4, whereas the impact of SRP, albedo, IR, and the navigation antenna thrust on the Galileo orbits is provided in Sect. 5. Section 6 contains a summary and conclusions.

2 Methodology

2.1 Box-wing model assumptions

In our analysis, we use the official metadata released by the European GNSS Agency. The metadata contains information about the satellite properties, i.e., the physical characteristics of both Galileo-IOV and FOC satellites such as mass, area, absorption, and reflectivity of the particular surfaces, as well as information about the attitude law and antenna parameters such as phase center offset and variations.

In this study, we use the box-wing model composed by Bury et al. (2019) which is formulated consistently with Rodriguez-Solano et al. (2012a) with the consideration of shadows resulting from both the Earth and Moon as a fraction of the eclipsed part of the Sun disk comprising antumbra and penumbra periods. The yaw-attitude corrections for low elevation of the Sun with respect to the orbital plane are not applied. The Galileo-FOC satellites are equipped with thermal radiators on \(-\)\(Z\), \(+Y\), \(-\)\(Y\), and \(-X\). During the yaw-steering, none of the Y-panels is illuminated by the Sun, nor is the −X panel which covers the clock. However, the delayed re-radiation which comes from the radiators may cause some effects. In this study, we consider only the immediate thermal re-radiation and neglect the remaining heating or cooling effects.

The optical properties of the satellite elements can change during their lifetime for the particular surfaces. The metadata contain two sets of values at the beginning and at the end of the satellite lifetime, which change for the satellite bus materials and are unchangeable in time for the solar panels. In our calculation, we assume the values which correspond to the beginning of the mission.

The box-wing model allows us to assess the characteristics of the accelerations resulting from both SRP and ERP which act on the Galileo satellites. The accelerations resulting from SRP have to be considered separately for the satellite bus and the solar panels. The satellite bus is covered by multilayer insulation for thermal protection. Hence, the accelerations \(a_{b}\) with immediate thermal re-radiation are described by the formula of Milani et al. (1987):

$$\begin{aligned}&a_{b} = -\frac{S_{c}}{c} \cdot \frac{A}{m} \cdot \cos \theta \cdot \nonumber \\&\quad \left[(\alpha +\delta )\cdot \left( e_{\odot }+\frac{2}{3}e_{n}\right) +2\rho \cdot \cos \theta \cdot e_{n}\right]\end{aligned}$$
(1)

In the case of the solar panels, we also assume the instantaneous re-radiation back into space; however, due to the fact that temperatures for the front and back sides are the same, and the thermal re-emission is mostly balanced, the accelerations \(a_{sp}\) resulting from SRP are described by formula:

$$\begin{aligned}&a_{sp} = -\frac{S_{c}}{c} \cdot \frac{A}{m} \cdot \cos \theta \cdot \nonumber \\&\quad \left[(\alpha +\delta )\cdot e_{\odot } + 2\left( \frac{\delta }{3}+\rho \cdot \cos \theta \right) \cdot e_{n}\right]\end{aligned}$$
(2)

In (1) and (2), \(S_{c}\) denotes the Solar constant (1367 W/\(\hbox {m}^{2}\)) rescaled by the \((\frac{AU}{r_{S}})^{2}\), where AU denotes the astronomical unit and \(r_{S}\) is an instantaneous distance of the satellite from the Sun, c is the speed of light, A stands for an area of a single flat surface element, m is the mass of the satellite. An angle between the unit vector of the surface normal \(e_{n}\) and the unit vector of the direction of the illuminating source \(e_{\odot }\) is described by \(\theta \). SRP results from the impulse transfer of the absorbed and emitted photons on the satellite’s surface illuminated by the Sun. Fractions \(\alpha \), \(\delta \), and \(\rho \) describe absorbed, diffusely reflected, and specularly reflected photons, respectively (with \(\alpha \) + \(\delta \) + \(\rho \) = 1, Milani et al. 1987).

Table 2 Characteristics of the Galileo orbit solutions
Full size table

In the case of the albedo and IR, the calculation of the perturbing acceleration is based on the mathematical formulae by Knocke et al. (1988) and is consistent with Rodriguez-Solano (2009); Rodriguez-Solano et al. (2012b) who performed a similar analysis on the impact of albedo on GPS Block IIA using Clouds and Earth’s Radiant Energy System data (CERES, Wielicki et al. 1996).

The Galileo metadata contain information only for the radiation in the visible spectrum. As a result, we use the same parameters for IR, which almost certainly are different to those describing the interactions in the visible spectrum.

Table 3 Characteristics of GPS Block IIF and Galileo satellites
Full size table
Fig. 2
figure 2

Accelerations acting on the particular GNSS satellites in direction D(left) and B(right) due to the direct SRP as a function of the angular height of the Sun above the orbital plane (\(\beta \)) and the argument of latitude of the satellite with respect to the argument of latitude of the Sun (\(\varDelta \)\(u\)). The scales for accelerations in the direction D and B are different due to significantly higher accelerations acting on satellites in the direction D. All acceleration are given in m/\(\hbox {s}^{2}\)

Full size image

2.2 Orbit determination strategies

We determine the Galileo orbit solution in several strategies using the modified version of Bernse GNSS Software 5.2 (Dach et al. 2015). We test analytical, empirical, and hybrid approaches.

For the absorption of SRP, the analytical approach (B) considers only the box-wing model. In that solution, we also apply both albedo and IR modeling, as well as the antenna thrust. In contrast, the empirical solution ‘E2’ does not consider the box-wing model for the SRP absorption and uses purely empirical parameters consistent with the ECOM2 model. In ‘E2,’ we consider ERP and antenna thrust. In order to concisely assess the impact of the neglected ERP and antenna thrust, we generate the solution ‘N2’ in which the ECOM2 parameters are responsible for the absorption of SRP and ERP. Hybrid strategies comprise the mixture of the analytical and empirical approaches for SRP absorption. We test three variants of the empirical parameters, that is, ‘H2’ in which we estimate a set of ECOM2 parameters, ‘H1’ in which we estimate five parameters of the ECOM1, and ‘H0’ in which we consider only constant accelerations in DYB directions. Additionally, we provide three solutions which test the impact of omission of albedo, IR, and antenna thrust, that is ‘EA2,’ ‘EI2,’ and ‘ET2,’ respectively. All the additional solution consider ECOM2 parameters without the box-wing models (Table 2).

GNSS observation processing is consistent with Zajdel et al. (2019) and Bury et al. (2019). We estimated a 1-day orbital arc integrating the orbits within 5-minute intervals. The solution is generated for the first 200 days of the year 2017 and is calculated based on the GPS and Galileo observations provided by approximately 100 globally distributed GNSS stations (Zajdel et al. 2019).

We analyze SLR residuals (Zhu et al. 1997; Urschl et al. 2007; Zajdel et al. 2017) for the particular orbit solutions as well as the orbit misclosures. We also check what values do empirical parameters assume when the box-wing model is applied, and finally, what is the impact of the particular forces, i.e., the SRP, albedo, IR, and antenna thrust, determined on the Galileo satellite positions.

3 Accelerations resulting from non-gravitational forces

3.1 Direct SRP

Based on the Galileo metadata, we calculate the perturbing accelerations which result from the direct SRP. The accelerations depend on the geometry of the SSE frame, i.e., the (\(\beta \)) and (\(\varDelta \)\(u\)). Table 3 contains the basic geometrical characteristics of the Galileo satellites compared to the GPS Block IIF spacecraft. We also plotted the accelerations resulting from the direct SRP which act on both GPS Block IIF and Galileo satellites in D and B directions as a function of \(\beta \) and \(\varDelta \)\(u\) (Fig. 2). We omitted the Y-direction due to the fact that, theoretically, no SRP accelerations should act on the Y-surface of the satellites. The calculated Y-accelerations for nominal yaw-steering did not exceed the magnitude of 10\(^{-11}\) m/s\(^2\) which translates into zero-level satellite position error. Figure 2 illustrates that the accelerations acting on the GPS satellites are significantly less dispersed than for the Galileo satellites.

Due to the fact that Galileo satellites are by a factor of more than two lighter than GPS satellites (see Table 3), they are more vulnerable to the direct SRP with the maximum absolute value of SRP accelerations at the level of 122 nm/s\(^2\) (excluding Galileo in eccentric orbits) as compared to 116 nm/s\(^2\) for the GPS satellites, in the direction D.

Fig. 3
figure 3

Amplitude spectra of the D—accelerations acting on GPS Block IIF (G06), Galileo-FOC (E08), Galileo-IOV (E12) and Galileo-FOC launched into highly eccentric orbit (E14) as a function of the Galileo satellite revolution period for different values of the \(\beta \)—angle

Full size image

The pattern of the plotted accelerations depends on the SSE geometry. The dispersion between different satellites in Fig. 2 comes from both the differences in Z-to-X surface area [corresponding to A from Eqs. (1) and (2)] ratio and the differences between the optical properties of surfaces. Galileo satellites are characterized by Z-to-X ratio at the level of 2.3:1 with a difference of 1.7 m\(^2\) as compared to GPS satellites which are almost cubic shaped spacecraft, i.e., the Z-to-X ratio for GPS spacecraft is at the level of 1:0.94 with a difference of 0.3 m\(^2\) (see Table 3).

Fig. 4
figure 4

Amplitude spectra of the B—accelerations acting on GPS Block IIF (G06), Galileo-FOC (E08), Galileo-IOV (E12) and Galileo-FOC launched into highly eccentric orbit (E14) as a function of the Galileo satellite revolution period for different values of the \(\beta \)—angle

Full size image

As a result, the maximum acceleration in D for the Galileo satellites occurs when both the \(\beta \) and \(\varDelta \)\(u\) are close to 0 because the Sun illuminates then the largest area −Z and whole solar panels, as well as during \(\beta \)\(\approx \)0 and \(-90^{\circ }~<~\varDelta \)u \(<~90^{\circ }\) when the Sun illuminates two surfaces, \(+X\) and \(-Z\). Moreover, the accelerations in the direction D are not purely symmetrical for the positive and negative \(\beta \) angles. This may stem from the fact that in this study we analyze the period of 200 days of 2017. The Earth orbits around the Sun on the orbit which is characterized by eccentricity e\(\approx \) 0.016. This results in the difference in the Sun–Earth distance at the level of 3.2% between aphelion and perihelion. As a result, the difference of the SRP acting on the Earth satellite between the two extreme points reaches the level of 6.4% because the SRP is rescaled by the change of the distance to the Sun raised to the power of two. Therefore, if our time series is only slightly longer than a half of the Galileo draconitic period, this explains the asymmetry of the plotted accelerations in the direction D.

The accelerations in the direction B are by two orders of magnitude lower than in the direction D for all the satellites. In the B direction for the Galileo spacecraft, the variations are significantly higher than for the GPS satellites, having the dominating nature of the once-per-revolution perturbation, which is confirmed by the spectral analysis (see Figs. 3 and 4). Periodic, nth-per-revolution \(n_\mathrm{ref}\) acceleration \(a_{Q}\) acting on the satellite in direction Q (e.g., in D, Y, B) can be expressed as follows:

$$\begin{aligned} a_{Q} = Q_{0}+ Q_{cn_\mathrm{ref}}\cos (n_\mathrm{ref}\cdot \varDelta u)+Q_{sn_\mathrm{ref}}\sin (n_\mathrm{ref}\cdot \varDelta u).\nonumber \\ \end{aligned}$$
(3)

When double integrating equation (3) and neglecting the integration constants which denote the mean offset and linear drift, one will obtain an amplitude \(\hbox {A}_{{m}}\) of the periodic error of the position of the satellite due to the particular periodic terms

$$\begin{aligned} A_{m}=\sqrt{Q_{cn_\mathrm{ref}}^{2}+Q_{sn_\mathrm{ref}}^{2}}\cdot \frac{1}{n_\mathrm{ref}^{2}}\cdot \left( \frac{T}{2\pi }\right) ^{2}. \end{aligned}$$
(4)
Table 4 Accelerations and satellite position error amplitudes due to the direct SRP for Galileo E08 FOC/E14 FOC ecc
Full size table

In Eq. 4, T denotes the satellite revolution period. Based on Eqs. (3) and (4), we calculated the periodic errors in the satellite positions and set them in Table 4 for Galileo-FOC E08 and E14. Moreover, we added information whether the particular terms are absorbed by the currently used ECOM2 model.

According to Table 4, the current set of the estimated ECOM2 parameters is suitable for the Galileo-FOC satellites. The highest periodic perturbations at the level of 371 and 146 mm, in B and D-directions, respectively, are absorbed by the ECOM2 model.

However, Galileo E14 and E18 are subject to the different perturbing forces due to the fact that they were launched into highly eccentric orbits. As a result, the non-considered once-per-revolution in D and twice-per revolution terms in B terms cause, at maximum, the errors at the level of 154 and 37 mm, respectively (see Table 4). The Galileo-IOV satellites behave similarly to the Galileo-FOC satellites due to the comparable mass and physical characteristics. GPS satellites are by the factor of 2 heavier than Galileo and have more favorable X-to-Z bus area ratio. As a result, the periodic perturbations are significantly smaller than for Galileo, i.e., the error caused by the once-per-revolution term in B is by a factor of 5.5 smaller for the GPS satellites. On the other hand, the spectral analysis indicates the same periodic perturbations; hence, the currently used ECOM2 is suitable for both Galileo in circular orbits and GPS satellites by means of absorbing the major first- and second-order SRP perturbations.

For the precise orbit determination with the accuracy at the 1–2 cm level, the standard ECOM2 is insufficient. The neglected thrice-per-revolution terms in the direction B cause perturbation at the level of 14 mm for the Galileo in circular orbits. Higher-order terms may interfere the orbit solution as well. However, estimation of the higher-order terms may destabilize the solution due to the increased number of parameters and the correlation between them (Rebischung et al. 2014; Bury et al. 2019). As a result, the application of the analytical a priori models such as the box-wing model may be beneficial for the solution.

3.2 Galileo vulnerability on the Y- and B-bias

Galileo satellites are designed to operate in the yaw steering mode apart from the eclipsing phase (Konrad et al. 2007; Montenbruck et al. 2015a). In the yaw steering mode, the satellite is constantly rotated about the Earth-pointing antenna boresight axis in the way that the solar panel rotation axis is perpendicular to both Earth and Sun directions (see Fig. 1).

As a result, we do not expect any periodic accelerations in the Y—direction, i.e., the calculated Y—accelerations for nominal yaw-steering did not exceed the magnitude of 10\(^{-11}\) m/s\(^2\) (see Sect. 3.1).

In the case of the GPS satellites, the Y-bias is well elaborated. Marquis and Krier (2000) explain that the Y-bias is caused by the internal heat re-radiated from the Y-surfaces; however, Fliegel et al. (1992) identified it with the misalignment of both solar panels and the bus w.r.t the direction of the Sun. Moreover, Kuang et al. (1996) calculated the misalignment of the solar panels w.r.t. nominal attitude for the GPS II satellites, obtaining deviation of the solar panels at the level of \(1^{\circ }\)–\(2^{\circ }\).

The Galileo-FOC satellites also suffer from Y-bias. This may be caused by the radiators mounted on each +/\(-Y\), \(-Z\), and \(-X\) panels. Moreover, the two radiators on Y panels are not symmetrical in terms of their area. This asymmetry, together with the non-instantaneously re-radiated heat from the radiators causes systematic effects. According to Sidorov et at. (2018), the effect caused by the radiator mounted on the \(-X\) panel can be mitigated by the estimation of \(D_{1S}\) term. However, the Y-bias has to be compensated by the estimation of an offset in the direction Y, that is by the \(Y_{0}\) term.

Fig. 5
figure 5

\(Y_{0}\) accelerations from the ECOM model acting on the particular Galileo satellites. Based on the CODE-MGEX products. The color of the line denotes a particular solution with different orbit modeling, whereas the gray planes indicate periods in which \(|\beta |\) < \(12.3^{\circ }\)

Full size image

We analyze the estimated Y-bias values based on operational CODE-MGEX products (Prange et al. 2017b). Figure 5 illustrates the Y—accelerations in the period 2014.0–2018.4 for the particular Galileo satellites. During the period, the orbit processing strategy has been changed several times. Before 2015 ECOM1 was used (Beutler et al. 1994; Springer et al. 1999). Since 2015, the ECOM2 proposed by Arnold et al. (2015) is used with the modification of the set of the estimated parameters. Before the latest upgrade of the orbit determination strategy, the Galileo-FOC satellites suffered from the modeling issues during the eclipsing period which is illustrated by the significant peaks in Fig 5. This pattern is mitigated with the latest strategy modification when CODE started using ERP, antenna thrust modeling, as well as activated the eclipse attitude law (Dach et al. 2018). However, no matter the strategy, for all the Galileo-FOC satellites a significant Y-bias is visible, whereas the IOV satellites are not affected. The systematic shift is at the level of \(-7.0\) \(\cdot \)10\(^{-10}\) m/s\(^{2}\) ± 0.6 \(\cdot \) 10\(^{-10}\) m/s\(^{2}\) for all the Galileo-FOC satellites; thus, the application of this systematic near-constant acceleration as an a priori value could possibly be considered because such a treatment may reduce the number of estimated parameters.

In contrast to Galileo-IOV satellites, the FOC spacecraft do not suffer from the \(\beta \) angle-dependent acceleration changes in the direction B. Figure 6 illustrates the estimated by CODE accelerations in the direction B, that is the constant coefficient \(B_{0}\), as a function of time. The highest magnitudes of the accelerations occur for the highest \(\beta \) angle and decrease when the altitude of the Sun descends with respect to the orbital plane. The mean value of the term \(B_{0}\) is at the level of −5.4 \(\cdot \)10\(^{-10}\) m/s\(^{2}\) for the Galileo-IOV E19. Lower values for the Galileo-IOV E12 come from the character of its orbital plane for which the maximum \(\beta \) angle values are lower than for E19, that is ± \(52^{\circ }\) and ± \(77^{\circ }\) for E12 and E19, respectively. According to Figure 6, the characteristics of the \(B_{0}\) term are insensitive to the SRP modeling strategy. The characteristic peaks during high \(\beta \) angle were visible when the classical ECOM1 was used as well as when ECOM2 was introduced. Moreover, the patter barely changed even when CODE started applying the ERP, antenna thrust, and the eclipse altitude modeling.

Fig. 6
figure 6

\(B_{0}\) accelerations from the ECOM model acting on the particular Galileo satellites based on the CODE-MGEX products. The color of the line denotes a particular solution with different orbit modeling, whereas the gray planes indicate periods in which \(|\beta |\) < \(12.3^{\circ }\)

Full size image
Fig. 7
figure 7

Accelerations acting on Galileo E08 in direction D(top) and B(bottom) due to the albedo (left) and IR (right) as a function of \(\beta \) and \(\varDelta \)\(u\). The scales for accelerations in the direction D and B are different due to significantly higher accelerations acting in the direction D

Full size image

Summarizing the sensitivity of the Galileo satellites to the accelerations in direction B and Y, the Galileo-IOV satellites suffer from the acceleration in the direction \(B_{0}\) which is dependent on the \(\beta \) angle. Galileo-FOC satellites, on the other hand, suffer from the constant Y-bias which is neither dependent on the \(\beta \) angle nor is sensitive to the difference in the orbit modeling, i.e., introducing modeling of albedo, IR, and antenna thrust. The only mitigated peaks in the Y-accelerations during the eclipsing period may result from the application of eclipsing attitude mode. We assume that the \(B_{0}\) accelerations for the IOV satellites, together with the Y-bias for the FOC satellites, may arise from the solar panel rotation lag, that is, the misalignment of solar panels with respect to the Sun direction (Guo et al. 2017) or the asymmetrical radiators distributed on different sides of the satellite bus.

3.3 Albedo and infrared radiation

The impact of the albedo reflectivity depends on the cloud coverage, the geographic location of the satellite, the seasonality, and the SSE geometry, whereas the IR acts on the satellite independently from the SSE geometry because the Earth emits the thermal radiation continuously.

Figure 7 illustrates the accelerations acting on the Galileo-FOC E08 due to albedo and IR. We discuss here only one Galileo-FOC satellite, due to the fact that for all the remaining spacecraft launched into near-circular orbit the accelerations are characterized by a similar nature. The accelerations are presented as a function of the \(\beta \) angle and \(\varDelta \)\(u\). Both albedo and IR cause the accelerations by two orders of magnitude lower than the direct SRP. In the direction D, the accelerations due to albedo occur when the satellite passes above the illuminated part of the Earth. The higher the \(\beta \) angle, the lower the accelerations. The highest value of the D-acceleration reaches 1.1 nm/s\(^{2}\) when \(\beta \)\(\approx \)\(0^{\circ }\).

In contrary to the albedo reflectivity (see Fig. 7), IR causes the negative perturbing accelerations in direction D for \(90^{\circ }\, {<}\, \varDelta \)\(u\,{<}\, 270^{\circ }\), due to the fact that thermal radiation is emitted also by the non-illuminated part of the Earth. As a result, the radiation act on the \(+Z\) panel of the satellite bus and the solar panels with a positive sign during \(\varDelta \)\(u\)\(\approx \)\(0^{\circ }\) and with a negative sign during \(\varDelta \)\(u\)\(\approx \)\(180^{\circ }\). The accelerations in the direction D reach the maximum values at the level of ±1.0 nm/s\(^{2}\) depending on the SSE geometry. As a result, the amplitude of the acceleration is higher for the infrared radiation as compared to the albedo reflectivity (see Fig. 7). The maximum value of acceleration due the infrared radiation in the B direction is at the level of 0.5 nm/s\(^{2}\).

Fig. 8
figure 8

Accelerations acting on Galileo E18 in the radial direction due to the summarized albedo and IR as a function of the satellite altitude (h) and \(\varDelta \)\(u\)

Full size image

According to Rodriguez-Solano (2009), the accelerations resulting from the ERP depend on the satellite altitude. The first two Galileo-FOC satellites orbit at the highly eccentric orbits; hence, their altitude above the Earth surface falls between 16,100 and 27,000 km. As a result, we can indicate the dependence of the accelerations resulting from both albedo and IR on the satellite altitude. Figure 8 illustrates accelerations resulting from both albedo and IR as a function of the satellite altitude (h) and \(\varDelta \)\(u\) for the radial component. The highest accelerations occur for \(\varDelta \)\(u\) close to \(0^{\circ }\), i.e., when the satellite passes above the surface of the Earth illuminated by the Sun. The remaining accelerations which occur when \(\varDelta \)\(u\) assumes values close to \(180^{\circ }\) result from IR. The higher the satellite altitude, the lower the magnitudes of the acceleration resulting from albedo.

Fig. 9
figure 9

Amplitude spectra of the D (top)—and B (bottom)–accelerations acting on Galileo-FOC E08 due to albedo reflectivity (left) and the infrared radiation (right) as a function of the Galileo satellite revolution period for different values of the \(\beta \)—angle

Full size image
Table 5 Accelerations and satellite position error amplitudes due to albedo and IR for Galileo-E08 FOC and Galileo-E14 FOC ecc
Full size table

We performed the spectral analysis of accelerations resulting from both albedo and IR (Fig. 9) in order to find characteristic periods by the example of the Galileo-FOC E08. Spectral analysis and the amplitudes of the periodic satellite position errors shown in Table 5 indicate which perturbations caused by albedo and IR are absorbed by the currently used ECOM2. The cm-level perturbations are visible for the once-per-revolution terms in both D and B directions. Since the odd terms are considered only in the B direction, the deficiencies in coping with the once-per-revolution acceleration causing errors at the level of 22 and 45 mm due to albedo and IR, respectively, have to be absorbed by the box-wing model. We also found the maximum accelerations in the radial component at the level of 1 nm/s\(^{2}\) for both albedo and IR which translates into a constant radial offset in the satellite position of 22 mm. In the case of the Galileo in eccentric orbits, the accelerations are slightly higher than for the Galileo in nominal orbits (see Table 5). The neglected in ECOM2 periodic acceleration in the direction D cause errors in the satellite positions at the level of 29 and 55 mm, due to albedo and IR, respectively.

Due to the fact that we use the same properties for the albedo and IR, the radial accelerations due to both forces are similar, with an error of 29 mm. In summary, both albedo and IR have to be considered from the a priori model due to the fact that ECOM2 cannot fully absorb the ERP-induced perturbations, yet it was designed purely for the SRP absorption.

4 Orbit determination results

We verify the internal and external accuracy of the solutions before analyzing the impact of the box-wing model application. The internal consistency of solutions is assessed based on the analysis of orbit misclosures and the external accuracy is evaluated based on the SLR residuals.

4.1 Evaluation of the internal accuracy

The consistency of the orbit solution is assessed based on the calculation of the satellite position differences at arc boundaries for each day decomposed into the radial, along-track, and cross-track directions. The orbit misclosures are calculated for the least stable parts of the orbital arc; hence, they comprise a reliable indicator for the evaluation of the consistency of the particular solutions.

Table 6 presents the results of the orbit misclosures calculated for all solutions. In general, all hybrid and empirical solutions are consistent at a similar level, with the exception of the purely analytical approach. The solution ‘B’ provides the least consistent orbits both in terms of the mean values and standard deviations (STDs) which shows that the purely analytical approach is insufficient in terms of the absorption of the direct SRP that has to be taken into account by a set of the empirical parameters. According to Fig. 10, the hybrid solution ‘H1’ is internally more consistent and accurate than the empirical ‘E2,’ despite a slightly higher mean value in the radial direction for ‘H1.’ Apart from that, the STD of the orbit misclosures for ‘H1’ is lower than for ‘E2,’ by 4, 2, and 1 mm for the radial, along-track, and cross-track component, respectively.

Table 6 Orbit misclosures for the Galileo-FOC satellites for the particular orbit solutions decomposed into the radial, along-track, and cross-track components. All the values are given in mm
Full size table
Fig. 10
figure 10

Galileo-FOC orbit misclosures decomposed into the radial, along-track, and cross-track directions for solutions ‘H1’ (magenta) and ‘E2’ (cyan). All the values are given in mm

Full size image

4.2 Evaluation of the external accuracy

The external accuracy of solutions is evaluated based on the SLR residual analysis. The SLR comprises a reliable, independent validation tool due to the fact that the SLR observations are not used in the precise orbit determination process as well as that the SLR uses wavelengths different to those employed in GNSS; that is optical instead of a microwave. Let us now scrutinize the processing results based on the analysis of SLR residuals as a function of the Sun elongation \(\varepsilon \) and \(\beta \) angles (see Fig. 1). The elongation angle is defined as:

$$\begin{aligned} \cos \varepsilon = \cos \beta \cos \varDelta u \end{aligned}$$
(5)

Figure 11 illustrates the SLR residuals for the FOC satellites for the particular solutions described in Sect. 2. The analysis is performed for Galileo satellites excluding the pair of the satellites launched into highly elliptic orbits. The red line in Fig. 11 indicates the linear dependence of SLR residual values on the elongation angle and the color of the dots indicates the value of the \(\beta \) angle. The statistics concerning the mean offset, the slope, STD, and STD values for the \(\beta \) angles below 12.3\(^{\circ }\), the eclipsing periods of Galileo satellites, are set in Table 7.

Fig. 11
figure 11

SLR residuals (vertical axis) as a function of the Sun elongation \(\varepsilon \) (horizontal axis) and \(\beta \) angles (color)

Full size image
Table 7 Statistics concerning the SLR validation for the Galileo-FOC satellites (excluding satellites on eccentric orbits) from particular solutions
Full size table

The solutions based on the empirical approach, that is ‘E2’ and ‘N2’ result in significant slopes at the level of −0.319 and −0.320 mm/\(^{\circ }\) for ‘E2’ and ‘N2,’ respectively. However, the solution ‘N2’ is affected by an offset higher by 27 mm than that from the solution ‘E2’ due to the omission of the antenna thrust and ERP.

The solutions based on the a priori box-wing model are characterized by a significantly lower slope. The slope for the solution ’B’ is at the level of 0.081 mm/\(^{\circ }\) with an offset at the level of −6.3 mm. However, the spread of the SLR residuals in higher than for the other hybrid solution. This indicates that a certain set of empirical parameters capable of absorbing the remaining accelerations has to be estimated. However, estimating seven of the ECOM2 parameters, we introduce a systematic effect due to the over parametrization of the solution which reveals in the slope at the level of −0.130 mm/\(^{\circ }\). In terms of the precise Galileo orbit determination, the solution based on the box-wing model with the reduced number of the estimated parameters seems to be the most suitable. Both the solutions ‘H0’ and ‘H1’ are characterized by almost flat trend of the SLR residuals, that is, 0.013 and 0.023 mm/\(^{\circ }\) for ‘H0’ and ‘H1,’ respectively. The STD of SLR residuals are below the level of 24 mm for both ‘H0’ and ‘H1’ solutions which is better by over 20% than in the empirical ECOM2 solutions. Moreover, the hybrid solutions are the most efficient during \(\beta \) angles below 12.3\(^{\circ }\). In the best case, the STD of SLR residuals for the eclipsing season is at the level of 22.1 mm for the solution ‘H1.’ Similar test was performed by Li et al. (2019) who applied the box-wing model with the classic ECOM1 model. However, Li et al. (2019) obtained the slope at the level of −0.17 mm/\(^{\circ }\). The only drawback of our box-wing solution is visible in the positive offset of the SLR residuals for all the hybrid and the analytical solutions for the FOC satellites. In the case of the IOV satellites (not shown here), the offset assumes the negative sign. This may stem from the fact, that we use the same properties for both albedo and IR, and we neglect the eclipsing attitude modeling.

To conclude, the box-wing model used together with estimating of a reduced set of the empirical parameters stabilizes the orbit solution and mitigates the systematic errors, especially during the eclipsing period.

5 Impact of non-gravitational forces on the Galileo orbits

5.1 Impact of the box-wing model application

Fig. 12
figure 12

Differences between SLR residuals to Galileo-IOV (top) and FOC (bottom) provided from solutions ‘H2’ and ‘E2’ as a function of the \(|\beta |\) and \(\varDelta \)\(u\). All values are given in mm

Full size image

The impact of the box-wing model application is assessed based on the differences between solutions ‘H2’ and ‘E2,’ because these solutions are different only in terms of the application of the box-wing model. Here, we analyze the differences between SLR residuals, satellite positions, as well as between the empirical parameters estimated in both solutions.

Figure 12 illustrates the differences between SLR residuals from solutions ‘E2’ and ‘H2’ as a function of \(\beta \) and \(\varDelta \)\(u\). In terms of coping with SRP especially crucial are the moments when the satellite eclipses. For the Galileo satellites on nominal orbits, the eclipse phase starts when the \(|\beta |\) is below \(12.3^{\circ }\), whereas for the Galileo on highly elliptic orbits the eclipsing phase initiates for \(|\beta |\) angles below 11.0–15.5\(^{\circ }\) depending on the orientation of the orbit perigee w.r.t the Sun. In Fig. 12, the highest differences are visible for \(|\beta |\) angles below \(20^{\circ }\). The current form of the ECOM2 model, which neglects the higher-order terms, is insufficient during the eclipsing season. However, ECOM2 was designed strictly for the yaw-steering mode, and for the \(|\beta |\) angle below \(4.1^{\circ }\) the Galileo switches to the ‘modified yaw-steering,’ which has to be taken into account by a special attitude modeling. The analytical box-wing model successfully absorbs most of the higher-order SRP for low \(\beta \) angles. The mean amplitude of the SLR residuals differences are at the level of 50 mm, and in the worst case, reaching 65 mm for the Galileo-FOC satellites (see Fig. 12—bottom).

However, the box-wing model on its own is insufficient for the complete absorption of the direct SRP, due to the fact that the external conditions, in which the satellite orbits, are challenging to model. This stems from the intensity variation of the photons which come from the Sun including solar winds as well as the continuous changes of the satellite attitude due to the yaw-steering and limitations in keeping the perfect satellite orientation. Moreover, the surface properties describing absorption, diffusion, and specular reflection are known with a limited accuracy and may also change after several months in orbit. As a result, it is almost impossible to control the perturbing accelerations when not estimating the constant acceleration in the direction D whose main task is to rescale the direct SRP effect.

Figure 13 presents absolute SLR residuals to the orbit solution based solely on the box-wing model for all Galileo-IOV as a function of \(|\beta |\) and \(\varDelta \)\(u\). In the SLR residuals to the Galileo-IOV orbits determined using the strategy ‘B’ we observe a systematic error for \(\varDelta \)\(u\) \(\approx ~90^{\circ }\) and \(270^{\circ }\) and \(\beta \) < \(60^{\circ }\). The reason for this pattern is unknown; however, it may be connected with a dependence of the accelerations in the direction \(B_{0}\) on the \(|\beta |\) angle. In solution ‘B,’ when the \(B_{0}\) term is not estimated, the errors in the satellite position increase above the absolute level of 50 mm, when \(\beta ~\approx ~90^{\circ }\) (see Fig. 13). In contrary to the IOV satellites, the residuals for the FOC satellite are free from the pattern observed in Fig. 13, because the FOC satellites are not affected by large B-bias.

Fig. 13
figure 13

SLR residuals to Galileo-IOV satellites derived from the solution ’B’ as a function of the \(|\beta |\) and \(\varDelta \)\(u\). The SLR residuals are given in mm

Full size image

Now, let us analyze the impact of the box-wing model application on the position of the satellites which are subjected to the third Kepler’s law of motion, that is Beutler (2004):

$$\begin{aligned} n^{2}a^{3} = GM \end{aligned}$$
(6)

where n denotes mean motion of the satellite, a stands for its orbit semi-major axis, and GM comprises the standard gravitational product. Hugentobler (2008) assumed that when changing the force model, that is of the GM parameter, by the application of the force in the radial direction \(R_{0}\) and not changing the mean motion of the satellite, one would obtain the reduction of the radial direction at the level of Sośnica (2014):

Fig. 14
figure 14

Differences in the radial component of satellite’s position from solutions ‘E2’ and ‘H2’ for the Galileo-FOC E08 as a function of \(|\beta |\) and \(\varDelta \)\(u\). All values are expressed in mm

Full size image
$$\begin{aligned} \varDelta r \approx -\frac{1}{3}\frac{a^{3}R_{0}}{GM} \end{aligned}$$
(7)

This can be observed in Fig. 14 illustrating differences of the satellite positions in the radial direction of Galileo-FOC E08. Owing to the application of the box-wing model, we can evaluate the impact of non-absorbed systematic effects which may arise from the higher-order terms of empirical SRP coefficients and other, non-considered thermal effects. The highest amplitudes of the position differences occur during the eclipsing season. The position differences are at the level from −80 to 65 mm, which is consistent with differences of SLR residuals from solutions ‘E2’ and ‘H2’ (see Fig. 12). Besides the effects related to the eclipsing period, Fig. 14 indicates an offset in the position differences at the level of 12 mm which also coincides with the mean value of the SLR residuals.

Eventually, we investigate what values do the empirical parameters assume when they are simultaneously estimated and the box-wing model is applied. According to Fig. 15, the box-wing model absorbs approximately 97% of the direct SRP acting in the direction D. Despite that, a residual acceleration at the level of −2.2 nm/s\(^{2}\) occurs. This acceleration has to be taken into account by the estimation of the SRP scaling parameter \(D_{0}\). The remaining 3% of the acceleration may result either from the limitations of the box-wing model, that is from the inadequate values describing properties of the satellite either from the inaccuracy of the Solar constant, or the deficiencies in the modeling of the satellite attitude. In the case of \(Y_{0}\) term, the Y-bias cannot be absorbed by the analytical models, therefore it assumes almost the same values in both ‘H2’ and ‘E2’ solutions.

Fig. 15
figure 15

Empirical parameters estimated in strategies ‘E2’ (left) and ‘H2’ (right) for the Galileo-FOC E09 expressed in nm/s\(^{2}\). The gray planes denote the \(|\beta |\) \(<~12.3^{\circ }\), i.e., eclipsing periods. Vertical scale for the \(D_{0}\) parameter (top) is different for the two solutions, due to the fact that accelerations in solution ‘E2’ are by two orders of magnitude higher than in solution ‘H2’

Full size image

The application of the box-wing model introduces a bias at the level of 0.82 nm/s\(^{2}\) in the direction B which is absorbed by the term \(B_{0}\). The periodic cosine terms \(D_{2C}\) and \(B_{1C}\) assume significantly lower values for the solution ‘H2’ which means that the respective accelerations are greatly absorbed by the box-wing model. The dependence on the \(\beta \) angle has been significantly diminished. However, the remaining accelerations estimated in the \(D_{2C}\) for the hybrid solution are characterized by a positive offset at the level of 0.72 nm/s\(^{2}\). The remaining accelerations in the \(B_{1C}\) term are characterized by a positive shift as well, but by the factor of two smaller than in \(D_{2C}\). The standard deviation of the accelerations is at the level of ±0.66 nm/s\(^{2}\). On the other hand, sine terms assume nearly the same values in both types of solutions. This means that the sine terms do not have the interpretation stemming from the box-wing model which assumes the instantaneous re-radiation of the heat coming from the SRP. As a result, the sine terms may absorb the accelerations which are shifted in phase with respect to the ‘cosine’ acceleration. This may come from to the absorption of accelerations arising from non-instantaneous heat re-radiation as well as for eventual solar panel lag. The standard deviation of the accelerations absorbed by terms \(D_{2S}\) and \(B_{1S}\) is at the level of 0.72 nm/s\(^{2}\) and 0.43 nm/s\(^{2}\), respectively.

Fig. 16
figure 16

The absolute median impact (orbit offset) due to the antenna thrust, albedo, and IR on the particular Galileo satellite types decomposed into the radial, along-track, and cross-track directions, respectively

Full size image

5.2 Impact of albedo, IR, and antenna thrust modeling

In order to evaluate the impact of albedo, IR, and navigation antenna thrust, we calculate additional orbit solutions, that is ‘EA2,’ ‘EI2,’ and ‘ET2,’ which omit the particular effects. Then, we calculated the differences of the Galileo satellite positions between the solution ‘E2,’ and the respective solutions.

The absolute values of the median differences (offsets) are illustrated in Fig. 16 for the Galileo-FOC and IOV satellites for the radial, along-track, and cross-track components. We also provide the STD statistics in Table 8 in order to describe the variability of the particular effects. According to Fig. 16, the radial component is the most affected by all the forces. The antenna thrust causes a constant radial acceleration due to the continuous broadcasting of the navigation signal. The effect for the Galileo-FOC satellites is higher than that for the Galileo-IOV satellites, i.e., 20.5 mm versus 13.5 mm, due to the fact that Galileo-FOC navigation antennas are more powerful, i.e., about 200 W for FOC satellites as compared to 155 W for Galileo-IOV. A similar effect was indicated by Steigenberger et al. (2017) who measured the Galileo-FOC antenna transmit power at the level of 265 W obtaining 26.5 mm offset which also translates into 10 mm offset effect per 100 W of the antenna transmit power. The accelerations resulting from the antenna thrust are stable during the whole satellite revolution for both FOC and IOV satellites due to the fact that these satellites orbit on nearly circular orbits. As described in Sect. 3.3, the amplitudes of the accelerations resulting from the IR are by the factor of 2 higher than those resulting from albedo. The results confirm these presumptions.

Table 8 Offset and STD of the median 3D impact of the particular forces acting on the particular types of the Galileo satellites
Full size table
Fig. 17
figure 17

Differences between SLR residuals to Galileo-IOV (top) and FOC (bottom) provided from solutions ‘N2’ and ‘E2’ as a function of the \(|\beta |\) and \(\varDelta \)\(u\). All values are given in mm

Full size image

The effect due to the albedo and IR omission in the radial direction reaches the level of 6.9 and 5.1 mm due to the albedo, and 16.1 and 14.2 mm due to the IR, for IOV and FOC satellites, respectively. The variability of the effect caused by IR is higher by almost the factor of two than albedo which results from the different characteristics of the two forces (see Sect. 3.3 and Table 8). Higher values of all effects, as well as their STDs for the Galileo-FOC on eccentric orbits, come from their orbit characteristic. The altitude on which the two eccentric satellites orbit falls between 16,100 and 27,000 km as compared to the nominal altitude of 23,000 km. As a result, the highest values of the effect caused by the albedo and IR occur when the satellites are the closest to the Earth. On the other hand, this has no impact on the maximum effect caused by the antenna thrust which remains at the same level as for the nominal Galileo-FOC satellites.

The combined impact of the albedo, IR, and antenna thrust was also checked by the analysis of differences between SLR residuals and the satellite positions from solutions ‘E2’ and ‘N2.’ Both solutions are based on the ECOM2 model and do not consider the a priori box-wing model.

Figure 17 illustrates the differences between SLR residuals from solutions ‘E2’ and ‘N2.’ In contrast to the effect caused solely by the SRP (Fig. 12), the difference of SLR residuals presented in Fig. 17 shows an offset at the level of approximately 30 mm due to the antenna thrust and IR which is consistent with Prange et al. (2017a). The impact of the IR increases when \(\varDelta \)\(u\) approaches \(0^{\circ }\) and \(180^{\circ }\) when the solar panel axis Y lies in the orbital plane, thus the irradiated area is the largest. The effect is only slightly lower for the Galileo-IOV satellites, however of the same characteristics as for the FOC satellites.

Table 9 The characteristic of the perturbing accelerations acting on the Galileo satellites
Full size table

6 Conclusion and discussion

In this study, we evaluate the impact of main non-gravitational forces acting on the Galileo satellites, that is the direct SRP, the Earth albedo, IR, and the navigation antenna thrust.

We found that the Galileo-IOV satellites suffer from accelerations in the direction B up to the level of −3.6 nm/s\(^{2}\) depending on \(\beta \) angle. On the other hand, the Galileo-FOC suffer from the Y-bias at the level of −0.7 nm/s\(^{2}\) which has not been seen for the Galileo-IOV generation. Differences between the accelerations for the two types of the Galileo satellites may originate from the fact that the FOC satellites are equipped with asymmetrical radiators on +/−Y surfaces as well as on the -Z and -X panels. Moreover, the accelerations may also occur due to solar panel rotation lag.

Based on the composed box-wing model, we assessed the accelerations resulting from SRP, albedo, and IR. The maximum values of the accelerations are listed in Table 9. In the case of the direct SRP, the neglected, odd periodic terms in the D direction have the largest impact on the satellites which orbit on the highly elliptic orbits (E14 and E18), i.e., the periodic errors caused by the once- and thrice-per-revolution accelerations are, at maximum, at the level of 154 and 23 mm, respectively. In direction B, for which the even terms are neglected, the highest errors occur for the twice-per-revolution accelerations and, at maximum, are equal to 37 mm. In the case of the albedo and IR, the biggest impact on the satellite positions have the once-per-revolution accelerations in the direction D which cannot be absorbed by the empirical ECOM2 model. The accelerations cause the maximum errors at the level of 22 and 45 mm for the Galileo-FOC E08 and Galileo-FOC E14 ecc., respectively. Most of the non-absorbed by the empirical models perturbations, together with higher-order perturbations, are to the great extent successfully absorbed by the composed box-wing model.

Additionally, we calculated that the acceleration due to the antenna thrust is at the level of −1.0 nm/s\(^{2}\) for the power of the antenna at the level of 200 W which translates into a radial satellite offset of about 20 mm.

In a series of different strategies, we estimate the precise Galileo orbits. Based on the analysis of the orbit misclosures, we found that all the solutions are consistent at a similar level, with the exception of the solution based solely on the analytical box-wing model. The solution based solely on the box-wing model (‘B’) is insufficient due to the unpredictable changes of the external conditions which ought to be absorbed by the estimation of the constant acceleration \(D_{0}\).

Based on the analysis of the empirical parameters from solutions ‘H2’ and ‘E2,’ we found that the box-wing model absorbs up to 97% of the SRP in the Sun-satellite direction. Moreover, we found a significant improvement of the orbit solution for the eclipsing period of the Galileo satellites when using the a priori box-wing model. This is also indicated by the differences between SLR residuals from solutions ‘H2’ and ‘E2’ whose STD are at the level of 50 mm, and reaches up to 65 mm for the Galileo-FOC satellites. The mitigation of the \(\beta \) angle-dependencies on the orbit solutions can be seen also when analyzing the empirical parameters from the solution ‘E2’ and ‘H2.’

The remaining 3% of SRP is absorbed greatly by the \(D_{0}\) parameter whose variability is at the level of 0.4 nm/s\(^{2}\). Moreover, the estimated values of the empirical parameters indicate that the sine terms \(D_{2S}\) and \(B_{1S}\) do not have an SRP-induced physical interpretation. The same applies to the accelerations in the Y-direction which are characterized by a systematic offset for the FOC satellites. Both sine terms, as well as B- and Y-biases, can be identified with the accelerations resulting from the thermal effect which occur due to non-instantaneous radiation from the radiators mounted on the satellites or the solar panel rotation lag.

Eventually, we evaluated the magnitude of the effect caused by the albedo, IR, and antenna thrust based on the differences between the SLR residuals from the solutions ‘E2’ and ‘N2’ and the satellites’ positions differences in the radial direction. The combined effect results in an offset at the level of 30 mm for all the Galileo satellites and increases to the level of 50–60 mm when \(\varDelta \)\(u\) reaches \(0^{\circ }\) and \(180^{\circ }\). The median impact of albedo, IR, and antenna thrust is at the level of 5, 14, and 20 mm, respectively, for the Galileo-FOC satellites.

The overall accuracy of the Galileo orbit solution was assessed based on the analysis of SLR residuals. The STD of SLR residuals for the hybrid solution ‘H0’ is at the level of 22.5 mm as compared to the standard ECOM2 solution for which we obtained STD of SLR residuals at the level of 28.7 mm for the Galileo-FOC satellites (excluding satellites on highly elliptic orbits). However, the improved standard deviation of the SLR residuals appears at the expense of the SLR residual offset which rises up to the level of 15.3 mm. This offset may arise either from a limitation of the box-wing model which includes inaccuracies of the optical parameters or deficiencies in the attitude modeling or usage of the same parameters for IR and albedo modeling. Despite all these issues, the box-wing model allows for increasing the accuracy of the Galileo 1-day orbit solutions to the level similar to that of GPS satellites from the long-term reprocessing (Sośnica et al. 2015).