Introduction

It has been shown that correct modeling of satellite attitude has a noticeable impact on Precise Point Positioning (PPP) applications [1]. In the case of GNSS satellites with large horizontal antenna offsets, improper attitude modeling affects both code and carrier phase observations, inducing decimeter level errors. The most important features and corrections that depend on correct attitude are [2]:

• Phase wind-up correction (PWC) • Antenna Phase Center Offset (PCO) • Solar Radiation Pressure (SRP) PWC describes the phase shift caused by the variation of antenna orientation, which can amount for errors of up to one narrow-lane cycle in ionosphere-free linear combination. PCOs refer to the differences between the satellite's position referred to its Center of Mass (COM), and the actual satellite antenna position where the GNSS signal originates. An incorrect yaw angle accounts for an error in computing the antenna position, which can in turn translate into a decimeter level loss in positioning accuracy [3]. Finally, SRP modeling is necessary for accurate orbits. In the case of GPS, inaccurate attitude modeling can increase errors caused by both solar and Earth radiation pressure [4].

Nominal attitude model

The attitude variations of the satellite during one revolution are governed by the need to orient the solar panels optimally towards the Sun, with the restraint of keeping the GNSS antenna oriented towards the center of the Earth. The nominal attitude model that best accomplishes this task is given in its most simplistic form by [5]:

Ψ 𝑛 = 𝐴𝑇 𝐴𝑁 2 − tan 𝛽 sin 𝜇(1)

where Ψ 𝑛 denotes the nominal yaw angle, 𝛽 (also called Sun beta angle) denotes the elevation of the Sun above the orbital plane, while 𝜇 is the angle between the satellite and the orbit midnight point (measured in the along-track direction). The orbit midnight is described as the point of the orbit farthest away from the Sun, and for 𝛽 angles between ±14.5°it passes through the Earth's shadow. In contrast, the orbit noon is defined as the point closest to the Sun. As the satellite passes the midnight, respectively noon points of the orbit, the nominal attitude model describes a rotation maneuver alongside the spacecraftfixed z-axis. Simply speaking, the satellite rotates while facing the Earth so that the solar panels keep tracking the Sun. These rotations occur solely as the satellite crosses the midnight and noon points, where 𝜇 becomes equal to 0°and respectively, 180°and are defined as midnight and noon maneuvers. But as 𝛽 approaches 0°, the nominal model presents a singularity during these maneuvers which cause the rotation rate of the satellite to increase to infinity, thus surpassing the physical capabilities of the satellite's Attitude Control System (ACS). In order to counteract this limitation, for low 𝛽 angles, the actual midnight and noon maneuvers deviate from the nominal model by employing lower, constant rotation rates. Each block of GPS satellites presents different attitude deviations, however in this work the focus will be on block IIF. Figure 1 showcases the satellite orbit with its noon and midnight points, 𝛽 and 𝜇 angles, alongside the spacecraft-fixed reference frame. In the considered frame, the axes are defined as:

• z-axis -direction of GPS transmitter antenna • y-axis -axis of rotation of solar panels • x-axis -vector perpendicular to the spacecraft-fixed YZ plane, pointing in the general direction of the Sun Additionally, the 𝛽 value measured on-board GPS IIF satellites contains a bias of −0.5°, introduced in order to counteract the effect of noise on the attitude control system when the satellite enters the shadow regime. Thus, the adapted model will have the final form [5]:

Ψ 𝑎 = 𝐴𝑇 𝐴𝑁 2 (︂ − tan 𝛽 sin 𝜇 )︂ + 𝐵(𝑏, 𝛽, 𝜇)(2)

𝐵(𝑏, 𝛽, 𝜇) = 𝐴𝑆𝐼𝑁

(︂ 0.0175 • 𝑏 sin 𝐸 )︂ , cos 𝐸 = cos 𝛽 • sin 𝜇(3)

where 𝑏 = −0.5°, and 𝐵 denotes the actual yaw bias introduced by 𝑏.

Reverse Kinematic Point Positioning

In what regards GPS satellites, there are a number of papers that address the attitude modeling of satellites pertaining to blocks IIA, IIR and IIF. Among these, the most exhaustive resources are found in [6], [7] and [8]. An exhaustive analysis of GPSS II/IIA and IIF maneuvers over 15.5 years can be found

x-offset (cm) y-offset (cm) z-offset (cm) Estimated 39.3 ± 1.9 -1.7 ± 1.3 127.4 ± 6.1 Manufacturer 39.4 0.0 109.3

Table 1

Estimated and nominal PCOs of GPS block IIF satellites, adopted from [8]. The estimated values were computed using the RKPP method.

in [9]. The standard method for attitude estimation of these satellites is the Reverse Kinematic Point Positioning (RKPP), introduced in [8], and briefly presented in [10] and [9]. The values in table 1 were adopted from [8]. The method is implemented by observing one satellite from multiple stations (receivers), and considering all relevant geodetic parameters fixed, including each of the station positions. These parameters are then used to estimate the satellite clock corrections and the antenna phase center positions epoch-wise. By knowing the horizontal antenna PCO, information about the satellite's yaw angle can be derived, provided that the satellite has a large enough horizontal antenna offset. The method is based on assimilating data from hundreds of receivers that have visibility over the spacecraft performing the maneuver. The RKPP method uses both phase and code ionospherefree combinations, that are additionally single-differenced between satellites to eliminate receiver clock offsets.

By using the RKPP method, the yaw angle of GPS IIF satellites was estimated and compared to the nominal model, and a series of deviations were found. Compared to the nominal model, the yaw angle of the real attitude of the satellites presents a linear drift as 𝛽 approaches 0°, alongside an occasional inversion of the maneuver rotation direction. This inversion was consistent with the −0.5°b ias introduced in the Sun angle. These characteristics differ between different blocks of GPS satellites, and this study focuses solely on the IIF satellites which present the same attitude model for the whole block. These findings were previously estimated using the RKPP method. The present study aims to implement a different method that estimates the same attitude of the GPS IIF satellites, based on phase residual signals. The residual signals are computed using a subnetwork of IGS stations, and use both code and phase observatios. The most notable difference between the present method and the RKPP method is that the number of stations used is of the order of a few dozens for the estimation of a single maneuver, as opposed to a few hundreds necessary for RKPP. The purpose of this study is to explain and test this method and show that the results obtained are similar to the RKPP method. The residual method tested in the next sections uses the stations azimuth and nadir angles to relate the residual to attitude deviations in the spacecraft-fixed XY plane, which in turn are used to estimate the yaw angle.

Residual attitude estimation method

In this section, a mathematical relation is derived in order to relate residual measurements to attitude deviations. The relation will then be adapted in the next sections to implement a set of linear equation systems for attitude estimation. In the following sections, vectors and matrices are represented by boldface letters, using lowercase notation for vectors, and uppercase notation for matrices, e.g. r ∈ R 𝑁 ×1 and A ∈ R 𝑁 ×𝑀 . Lowercase italic letters denote scalar numbers, e.g. 𝑟 ∈ R. Subscript indices denote reference to a certain receiver, e.g. 𝜌 1 refers to the residual signal computed for station 1, while superscript indices denote reference to a satellite, e.g. 𝜏 𝑠 refers to the time correction of satellite 𝑠. It is also assumed that the italic expression of a (boldfaced) vector denotes its absolute value, e.g. 𝑟 = |r|.

Modeling the relation between attitude and residuals

Phase residuals necessary for yaw estimation are generated using the Bernese GNSS Software V5.2 [11], while the raw data used to obtain the residuals was downloaded from the International GNSS Service (IGS) archive in the form of RINEX files with measurements taken once every 30 seconds. Initially, a sub-network of IGS stations needs to identified which has visibility over the maneuvering satellite. Utilizing a PPP procedure along with orbit and clock products obtained from the Center for Orbit Determination in Europe (CODE) [12], [13], the coordinates of each station are computed. Along with these coordinates, the receiver clock corrections and troposphere parameters are computed as well. After this first step, the precise coordinates of each of the desired stations alongside the aforementioned parameters will be considered known and fixed. This step excludes the maneuvering satellite from the calculations, since its attitude deviations would otherwise impact the initial precise estimation of the geodetic parameters.

Once the station coordinates are established, the focus shifts to the maneuvering satellite, and the previously estimated parameters are now considered fixed. In this next step, phase ambiguities and satellite clock corrections are estimated, considering all station-specific parameters fixed (station coordinates, receiver clock corrections, troposphere parameters), as they were estimated based on a PPP procedure during the first step. When computing the residual signals, solely the satellite antenna z-offset is accounted for, leaving the horizontal offset alone. As a result, the obtained phase residuals contain essential information regarding the satellite's antenna offset in the body-fixed x and y directions, as well as the yaw angle evolution during the maneuver. This enables accurate estimation and analysis. The residual can be expressed as:

𝜌 = 𝑑 ′ − 𝑑 = |r + p − r r | − |r − r r |(4)

where we define the following variables: Since the nominal satellite position is computed with regard to its COM, and the real origin point of the GNSS signal is in the center of the transmitter antenna, the real satellite's position has to be corrected by p. It can also be noted that p expressed in figure 2 is in the XY plane of the spacecraft-fixed reference frame. This is an approximation, since the z-axis offset amounts to a very small effect on the residuals, as absorbed by the estimated phase ambiguities and satellite clocks. Therefore, from now on p will be used to denote the projection of the antenna offset vector in the orbit-fixed XY plane. Finally, the residual signal becomes:

𝜌 = |r + p − r r | − |r − r r | = r − r r |r − r r | • p = e • p = −e ′ • p(5)

where e = r−rr |r−rr| denotes the unit vector pointing from the receiver to the satellite, while e ′ = −e.

Furthermore, the residuals can be expressed in terms of the yaw, azimuth and nadir angles. The yaw angle is denoted by Ψ -the angle between the velocity vector v and the x-axis of the spacecraft-fixed reference frame. The nadir angle denoted by 𝜂 is the angle under which the satellite observes the receiver with respect to the vertical direction. The azimuth angle is denoted by 𝛼. In figure 2 it is expressed as the angle between the velocity vector v and the projection of the direction vector e ′ onto the XY plane, denoted as e ′ 𝑥𝑦 . The nadir, respectively azimuth angles are the polar angles of the station position expressed in the spacecraft-fixed reference frame, and their derivations are introduced in Appendix A. Then, the displacement (offset) vector p and the direction vector e ′ can be expressed:

p = 𝑥 0 ⎛ ⎝ cos Ψ sin Ψ 0 ⎞ ⎠ e ′ = ⎛ ⎝ cos 𝛼 sin 𝜂 sin 𝛼 sin 𝜂 cos 𝜂 ⎞ ⎠(6)

where 𝑥 0 denotes the x-axis nominal antenna offset, as given in table 1. By inverting the sign convention in (5) for simplicity, the residual signal can be expressed as:

𝜌 = p • e ′ = 𝑥 0 (cos Ψ cos 𝛼 + sin Ψ sin 𝛼) sin 𝜂 = 𝑥 0 cos Ψ • sin 𝜂 cos 𝛼 + 𝑥 0 sin Ψ • sin 𝜂 sin 𝛼 (7)

Once the residuals are known, one can use an Least Squares Estimation (LSE) procedure to estimate the unknown terms 𝑥 ˚and 𝑦 ˚, and compute the yaw angle Ψ using the terms derived from (7):

𝑥 ˚= 𝑥 0 • cos Ψ , 𝑦 ˚= 𝑥 0 • sin Ψ , Ψ = 𝐴𝑇 𝐴𝑁 2 (𝑦 ˚, 𝑥 ˚)(8)

Estimating attitude using residuals

In this section, a set of estimation techniques will be presented, which solve different equation systems for the yaw angle using phase residuals. First, the unconstrained estimation method will be explained. Then, a series of models that solve the aforementioned equations will be showcased, introducing a linear constraint and estimations of satellite clock biases. In the following section, the equation systems are time-dependent with variables changing at each epoch. For simplicity, the time dependence will not be explicitly included in the notation.

Estimation of epoch-wise yaw angle

If the x and y biases introduced by the yaw angle are denoted by 𝑥 ˚(𝑡) and 𝑦 ˚(𝑡), the phase residual sample at epoch 𝑡 can be expressed as 𝜌(𝑡) = 𝜌(𝑥 ˚(𝑡), 𝑦 ˚(𝑡)), and the problem can be solved for these two variables accordingly. If a large enough number of observations are available, then a linear system of equations can be constructed for each time epoch. Then, the yaw angle can be estimated by relating the residual value to the biases introduced in the spacecraft-fixed XY plane. The epoch-wise equation systems can be solved using an LSE procedure.

The system is defined as 𝜌 = A • p, and the observations vector is expressed as:

𝜌 = [︀ 𝜌 𝑠 1 (𝑡) 𝜌 𝑠 2 (𝑡) . . . 𝜌 𝑠 𝐾 (𝑡) ]︀ 𝑇(9)

The observation vector 𝜌 is unique for each time epoch, and it contains every residual sample 𝜌 𝑠 𝑘 (𝑡) generated by each station observing satellite 𝑠 at time epoch 𝑡. The unknowns vector is denoted as p and contains the two components of the bias introduced by the unknown yaw angle Ψ(𝑡) expressed in the XY plane. The p vector is also time-dependent:

p = [︀ 𝑥 ˚(𝑡) 𝑦 ˚(𝑡) ]︀ 𝑇(10)

For each epoch, there should be at least two observations available, meaning that the satellite should be visible from at least two stations at all times. The unknown biases 𝑥 ˚(𝑡) and 𝑦 ˚(𝑡) will be solved for and then used to compute the yaw angle Ψ(𝑡) per epoch, according to (8). The term 𝜌 𝑠 𝑘 (𝑡) will be denoted as 𝜌 𝑘 , for simplicity. Finally, the grand design matrix is defined as:

A = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 𝜕𝜌1 𝜕𝑥 ˚𝜕𝜌1 𝜕𝑦 ∂𝜌2 𝜕𝑥 ˚𝜕𝜌2 𝜕𝑦 .

. . . . .

𝜕𝜌 𝐾 𝜕𝑥 ˚𝜕𝜌 𝐾 𝜕𝑦 ˚⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎝ sin 𝜂 1 cos 𝛼 1 sin 𝜂 1 sin 𝛼 1 sin 𝜂 2 cos 𝛼 2 sin 𝜂 2 sin 𝛼 2 . . . . . . sin 𝜂 𝐾 cos 𝛼 𝐾 sin 𝜂 1 sin 𝛼 𝐾 ⎞ ⎟ ⎟ ⎟ ⎠(11)

where 𝜂 𝑘 = 𝜂 𝑠 𝑘 (𝑡) and 𝛼 𝑘 = 𝛼 𝑠 𝑘 (𝑡) denote the nadir and azimuth angles under which satellite 𝑠 observes station 𝑘 at epoch 𝑡. Finally, the values of 𝑥 ˚and 𝑦 ˚at epoch 𝑡 are obtained by computing the pseudo-inverse: p =

(︀ A 𝑇 A )︀ −1 A 𝑇 • 𝜌

Constraining antenna offset

The previous section described a simple model for estimating the yaw angle using only acquired observations and nadir, respectively, azimuth angles. However, due to measurement noise, this estimation is not always precise. It is possible to improve this by adding a constraint [14], and forcing the absolute value of the horizontal antenna offset to 𝑥 0 :

min 𝑝 ‖Ap − 𝜌‖ 2 2 s.t. Cp = d(12)

where Cp = d is a linear constraint. Then the system of equations introduced above can be solved using the Karush-Kuhn-Tucker equations [15]. The constraint is introduced as:

𝑥 ˚2 + 𝑦 ˚2 = 𝑥 2 0 • (︀ cos 2 Ψ + sin 2 Ψ )︀ = 𝑥 2 0 (13)

This relation applies to the measured x and y-axis biases and constrains them on the circle described by the antenna rotation. However, this system requires an iterative approach since the constraint is non-linear and uses previously computed solutions for 𝑥 ˚and 𝑦 ˚. The initial values can be acquired by solving the unconstrained problem. However, due to the constraint being non-linear, a reformulation is necessary, which is introduced in Appendix B. The final form of the optimization problem with a linear constraint is derived as:

min 𝑝 ‖Ap 𝑗+1 − 𝜌‖ 2 2 s.t. [︀ 𝑥 ˚𝑗 𝑦 ˚𝑗]︀ • p 𝑗+1 = 1 2 (︀ 𝑥 2 0 + 𝑥 ˚2 𝑗 + 𝑦 ˚2 𝑗 )︀(14)

where

p 𝑗+1 = [︀ 𝑥 ˚𝑗+1 𝑦 ˚𝑗+1

]︀ 𝑇 denotes the solution computed for time epoch 𝑡 at iteration 𝑗 + 1. The term 𝑥 0 ≈ 39.4 cm denotes the antenna x-axis offset, i.e., the radius of the antenna rotation circle, as presented in table 1.

Estimation of satellite clock errors

Furthermore, it can be taken into account that some of the errors caused by poorly modeled attitude were absorbed by satellite clock corrections into the bias term 𝛿𝜏 𝑠 ˚, and are no longer included in the residuals. This was also observed and studied in [16]. In this case, the residual equation ( 7) can be altered to contain this bias and include it in the estimation:

𝜌 = 𝑥 0 cos Ψ • sin 𝜂 cos 𝛼 + 𝑥 0 sin Ψ • sin 𝜂 sin 𝛼 + 𝑐𝛿 ˚𝜏 𝑠(15)

After applying this change, the grand design matrix and unknowns vector change to: while the observations vector remains as in (9). For this model, at least three stations should observe the satellite at all times during a maneuver. Since the satellite clock offset bias is satellite-related and only one satellite is observed, there will only be one value per epoch to estimate for all stations.

A = ⎛ ⎜ ⎜ ⎜ ⎝ sin 𝜂 1 cos 𝛼 1 sin 𝜂 1 sin 𝛼 1 1 sin 𝜂 2 cos 𝛼 2 sin 𝜂 2 sin 𝛼 2 1 . . . . . . . . . sin 𝜂 𝐾 cos 𝛼 𝐾 sin 𝜂 1 sin 𝛼 𝐾 1 ⎞ ⎟ ⎟ ⎟ ⎠ , p 𝑗+1 = ⎛ ⎝ 𝑥 ˚𝑗+1 𝑦 ˚𝑗+1 𝑐𝛿 ˚𝜏 𝑠 𝑗+1 ⎞ ⎠(16)

Finally, the constraint is reformulated as in (12), with

C = (︀ 𝑥 ˚𝑗 𝑦 ˚𝑗 0 )︀ and 𝑑 remains a scalar, 𝑑 = [︀ 𝑥 2 0 + 𝑥 ˚2 𝑗 + 𝑦 ˚2 𝑗 ]︀ /2.

Then the optimization problem is formulated using the updated grand design matrix A and unknowns vector p 𝑗+1 from ( 16).

Results

In this section, a set of 70 maneuvers will be analyzed, spanning over two eclipse seasons between day 350 of the year 2020 and day 65 of 2021. In this time frame, four GPS IIF satellites were included in the analysis, namely PRN 25, 26, 27, and 32. These satellites perform 4 maneuvers per day: two noon maneuvers (abbreviated by N1 and N2) and two midnight maneuvers (abbreviated by M1 and M2). Furthermore, an implementation of the residual estimation method is made available under the corresponding author's Github page. The same code was used to generate the following figures and results.

In order to compute the following results, residual data was generated for each of the four satellites during their maneuvering times. As expressed Section 2, the phase residuals used in this analysis were generated using code and phase observations acquired from the IGS data archive, alongside satellite orbit products generated by CODE. The clock corrections and maneuver times were computed based on CODE orbit products, and only phase residuals were considered for the final attitude estimation.

Each maneuver presented in the following sections was observed by around 30 IGS stations, a number lower than the standard for RKPP. Furthermore, the residual data was screened in two different ways. The first data set contains the raw residuals from all stations that had visibility of the respective satellite during their performed maneuver. The second set of residual data was additionally filtered using a low-pass Butterworth filter of order 8 with a cutoff frequency of (1/4)𝑓 𝑠 = 1/120 Hz, adapted to the sampling frequency of the residual signals, at two samples per minute. This helped to further minimise effects of noise and outliers, and isolate the slow variations of the residual signals. In figure 3 two example sets of residual signals are shown. The data was captured during two noon maneuvers, and each line represents residuals captured by one IGS station observing the maneuvering satellite. The maneuvers start and ending times can be identified by the bump present in the residual signals. The residuals presents a noticeable increase during ongoing maneuvers, caused by attitude deviations from the nominal model, as observed previously in [8] and [16].

Linear drift

The linear drift of the yaw angle during noon and midnight maneuvers was observed in multiple studies. Basically, for low sun angles with an absolute value of |𝛽| ≤ 4.5 °, the nominal rotation rate becomes too fast for the ACS, and the satellite starts the maneuver by rotating at a constant rate until the maneuver ends and nominal attitude is resumed. In [8] this rotation rate is observed to have a nearly constant value of ≈ 0.06 °/𝑠 for midnight maneuvers, while in [7] a rate of ≈ 0.09 °/𝑠 − 0.13 °/𝑠 is observed for noon maneuvers. Similar findings were computed in this study, and are described in the following paragraphs.

Figure 6 illustrates the estimated rotation rates for the entire set of 70 maneuvers, and it becomes clear that the rotation rate deviates from the nominal model: as 𝛽 approaches zero, the nominal rotation rates present very high values, while the real, estimated rates remain low. In figure 7, noon and midnight maneuvers are separated, making it clear that noon turns present a slightly higher rotation rate than midnight turns. A set of examples of noon maneuvers can be observed in figure 4. In red, the nominal model given by ( 1) is shown, alongside the estimated yaw angle (in blue). Every previous study found that during noon turns, GPS IIF satellites start the maneuver at roughly the same time as indicated by the nominal model, however the rotation rate is not fast enough to catch up with that of the nominal model. The rotation rate therefore becomes constant, at an estimated value of ≈ 0.06 °/𝑠, meaning also that it takes a considerably longer time for the maneuver to be performed. After the turn is completed, sometimes a slight post-turn deviation can also be observed for the noon maneuvers, which is corrected as the satellite resumes its nominal attitude regime.

Similarly, the midnight turns from figure 5 also present similar deviations: the maneuver again has a constant rotation rate at a value of ≈ 0.088 °/𝑠 − 0.1 °/𝑠, slightly higher than the one observed during noon turns. Moreover, the midnight maneuvers start at a different time than the nominal model indicates: as can be observed from the figures, the satellite starts rotating directly after shadow entry, instead of performing the maneuver when reaching the middle region (the midnight point). The maneuver ends also visibly later than the nominal model indicates, just as the satellite exits the shadow region. All previously mentioned deviations were observed in numerous studies beforehand.

Table 2 presents the summarized rotation rates estimated in this study over the whole set of 70 maneuvers for the four satellites. The first line depicts results obtained based on the raw residual data, while the second line shows results obtained based on data filtered with the low pass Butterworth filter that helped to get rid of high-frequency variations. The estimated rotation rates were separated into percentiles, to further isolate the influence of outliers, and average and median values were computed for each sub-set. According to the data summarized in the table, the rotation rates for midnight maneuvers

Yaw bias effects

Another feature that is not always properly modeled when estimating yaw values is the effect of the bias applied to the observed Sun angle. It's effect on attitude maneuvers was observed in [16] and in [5].

As described in (3) the Sun angle value is biased by −0.5°in order to control the steering of the satellite during noon and midnight turns when 𝛽 = 0°. This causes an inversion of the sign of the maneuver performed under certain low 𝛽 values. In the presence of a yaw bias the flip of the direction of satellite rotation during yaw maneuver no longer appears at 𝛽 = 0°, i.e., when the Sun crosses the orbital plane, but for a negative 𝛽 angle. This thus leads to an inversion of the sign of the rotation rate with respect to the nominal rate for low 𝛽 values.

In the present data set, six noon maneuvers were observed at −1°≤ 𝛽 ≤ 0°. It was expected that out of these maneuvers, those performed at −0.5°≤ 𝛽 ≤ 0°would be inverted with respect to the nominal rotation direction [16]. In addition to that, it was expected that those maneuvers performed at 𝛽 between −0.9°and roughly −0.5°might occasionally exhibit an inverted rate, as expressed in [17] and [7]. This is caused by noisy estimations of 𝛽 on-board the satellite. All noon maneuvers performed with −1°≤ 𝛽 ≤ 0°are logged in table 3. The "Sign" column indicates whether the maneuver is inverted (-) or if it follows the direction of the nominal model (+). Within this data set there are four inverted maneuvers, out of which three are performed under a mean angle of −0.5°≤ 𝛽 ≤ 0°. One of the inverted maneuvers however does not correspond to the expectation, since it is performed under an angle of 𝛽 = −0.9°. It is possible that this is a satellite-specific behaviour for PRN 25, but in any case the present data set contains only a small number of inverted maneuvers, thus making it impossible to analyze the behaviour of a specific satellite individually.

As a last step, all the midnight maneuvers performed with −1°≤ 𝛽 ≤ 0°are also summarized in table . In [17] it was stated that the inversion caused by the −0.5°bias affects exclusively noon maneuvers, meaning that midnight maneuvers do not present inverted rotation rates. At the same time, [16] notes that midnight maneuevrs are sometimes inverted as well. In this study, one of the analyzed midnight maneuvers performed a rotation in the inverted direction, contrary to the analysis presented in [17] and according to the findings presented in [16]. However, none of the other midnight maneuvers were inverted. This midnight maneuver can be seen in figure 8, and was performed by PRN 26 at a mean Sun angle of 𝛽 = −0.05°.

Corrected residuals

After analyzing the effects of linear drift on the rotation rates and effects of the Sun angle bias on satellite attitude, a further analysis was made using the estimated yaw angles to correct residual signals.

As observed in figure 3, the rotation of the antenna around the z-axis during a midnight or noon turn generates an increase in residual signals which then fades away as the yaw angle approaches nominal values again. According to the model introduced in ( 15), the residual signal can be corrected not only for yaw values, but also for biases induced by clock deviations.

Clock corrections

Finally, the effect of estimating the biases on clock corrections was also analyzed. It is assumed that wrong modeling of the antenna phase center position can partially be absorbed by the estimated satellite clock corrections, thus causing abrupt variations during maneuvers. Figure 10 showcases two sets of satellite clock corrections computed during noon maneuvers for PRN 25. Close to the middle of the figure, as the maneuver is performed, a very abrupt dip in the clock correction can be observed. The original clock products are obtained using CODE products [18]. The estimated clock corrections are illustrated in figure 10 and are conclusive with findings presented in [16]: applying the estimated clock biases on the original clock corrections has the effect of reducing the abrupt change at the beginning and end of the performed maneuver. This results in smoother curves of the satellite clock corrections, and lower magnitudes of the spikes present during both maneuvers.

Conclusion

In the previous sections, a new approach to estimate the attitude of GPS IIF satellites was presented. The algorithms were applied to a set of residual data spanning over approximately 70 maneuvers, including four satellites in the analysis. The results obtained using the residual estimation method were analyzed with regard to linear drift, yaw bias effects on rate inversion, as well as effects of corrected yaw modeling on residuals and satellite clock corrections.

All aforementioned aspects were found to coincide with findings from all previous studies, indicating that the presented attitude estimation method is valid and yields good results. However, it was also visible that the method is rather sensitive to measurement noise present in the residual signals. This caused the estimated yaw angles to contain rather large outliers, which, in turn, affected the rotation rate estimations and occasionally the satellite clock and residual corrections presented in the final sections. This indicates that there is a need to develop a more advanced filtering method to eliminate outliers and measurement noise from residual signals before performing the yaw angle estimation.

Appendix A: Expressing azimuth and nadir angles

This section explains how the unknown angles for azimuth, 𝛼, and nadir, 𝜂, are derived by knowing the satellite and receiver positions. As expressed beforehand, these angles are related to the observed station position in the spacecraft-fixed reference frame. These angles can be expressed as functions of the satellite and receiver position, as well as the satellite's velocity vector. Next, the angular momentum vector h and its length are defined. The angular momentum vector h is perpendicular to the instantaneous orbital plane defined by v and r, as observed in figure 2. Then this can be used to compute the projection of the receiver's position vector r r , onto the plane normal to the nadir direction, r r⊥ . If a circular orbit is assumed, as is approximately the case for GPS satellites, the following equations hold:

h = r × v , ℎ = 𝑟 • 𝑣 , r r⊥ = 𝑟 𝑟 cos 𝛼 • v 𝑣 + 𝑟 𝑟 sin 𝛼 • h ℎ(17)r r⊥ = r r − (︁ r r • r 𝑟 2 )︁ • r(18)

where 𝑣 and ℎ define the lengths of the velocity, respectively angular momentum vectors. Per definition, the angular momentum vector h is also located in the spacecraft-fixed XY plane. In the case of circular orbits, v ⊥ r, meaning that v is also in the XY plane. Since GPS orbits are approximately circular, this allows for the yaw angle Ψ to be measured as the angle between the spacecraft-fixed x-axis and the velocity vector v. Furthermore, the product between r r⊥ and v, respectively h, can be expressed:

r r⊥ • v = 𝑟 𝑟 ⊥ • 𝑣 cos 𝛼 (19)r r⊥ • h = 𝑟 𝑟 ⊥ • ℎ sin 𝛼 = 𝑟 𝑟 ⊥ • 𝑟𝑣 sin 𝛼(20)

This then leads to the formula for computing the azimuth angle, 𝛼, and the nadir angle, 𝜂:

𝛼 = 𝐴𝑇 𝐴𝑁 2 (︂ r r⊥ • h 𝑟 , r r⊥ • v )︂ , cos 𝜂 = − (r r − r) • r |r r − r| • 𝑟 (21)

Appendix B: Derivation of linear constraint of antenna offset

The constraint of the antenna movement on a circle with constant radius 𝑥 0 was previously described in (13). If an iterative approach is considered, the estimated parameters 𝑥 ˚and 𝑦 ˚for iteration 𝑗 + 1 can be approximated linearly as: where Δ𝑥 ˚and respectively Δ𝑦 ˚denote the unknown corrections applied at iteration 𝑗 + 1 to the solutions 𝑥 ˚𝑗 and 𝑦 ˚𝑗. Then the constraint can be reformulated, by replacing the terms in (13) where 𝑥 ˚𝑗 and 𝑦 ˚𝑗 denote the known variables, estimated for the previous iteration. The final form of the constraint will then become:

(︀ 𝑥 ˚𝑗 𝑦 ˚𝑗)︀ • (︂ 𝑥 ˚𝑗+1 𝑦 ˚𝑗+1 )︂ = 12(︀ 𝑥 2 0 + 𝑥 ˚2 𝑗 + 𝑦 ˚2 𝑗 )︀(28)

with 𝑥 0 ≈ 39.4 cm denoting the radius of the antenna rotation circle, as presented in table 1.