The scope of this work is to present a new approach for estimating the attitude of GNSS satellites with horizontal ofset antennas. Previously, studies have shown that the attitude deviation from nominal models can successfully be measured by Reverse Kinematic Point Positioning (RKPP) approaches. Here, a new approach is studied which uses ionosphere-free observation residuals to estimate the yaw angle of these satellites. This approach is analyzed and applied to a two months set of data generated by observing GPS block IIF satellites. The method is derived and tested against previous works obtained using the RKPP method and the results coincide with findings of previous studies. Features such as the linear drift of yaw maneuvers during eclipse seasons, observed inversion of maneuvers performed under certain low Sun angles, as well as drifting of rotation rates are measured and their values are shown to be consistent with previous findings.
It has been shown that correct modeling of satellite attitude has a noticeable impact on Precise Point
Positioning (PPP) applications [
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
diferences 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 [
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 [
y
x z µ
ß Midnight
Noon
Additionally, the value measured on-board GPS IIF satellites contains a bias of − 0.5°, introduced
in order to counteract the efect of noise on the attitude control system when the satellite enters the
shadow regime. Thus, the adapted model will have the final form [
+ (, , ) ︂( 0.0175 · )︂
sin (, , ) = cos = cos · sin (2) (3) where = − 0.5°, and denotes the actual yaw bias introduced by .
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
[
Manufacturer
in [
The values in table 1 were adopted from [
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° bias introduced in the Sun angle. These characteristics difer between diferent 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 diferent 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 diference 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.
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|.
Phase residuals necessary for yaw estimation are generated using the Bernese GNSS Software V5.2 [
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-ofset is accounted for, leaving the horizontal ofset alone. As a result, the obtained phase residuals contain essential information regarding the satellite’s antenna ofset 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 − rr| − | r − rr| (4) where we define the following variables: • r - position vector of the satellite, related to its COM projection in the body-fixed XY plane • rr - the station’s position vector • p - antenna ofset vector in body-fixed x and y directions, as in figure 2 • - nominal distance between station and satellite, = |r − rr| • ′ - actual distance between station and the satellite’s antenna phase center, ′ = |r + p − rr| 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 ofset amounts to a very small efect 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 ofset vector in the orbit-fixed XY plane. Finally, the residual signal becomes: = |r + p − rr| − | r − rr| = |rr −− rrrr| · p = e · p = − e′ · p (5) where e = |rr−− rrrr| 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 (ofset) vector p and the direction vector e′ can be expressed: ⎛cos Ψ⎞ p = 0 ⎝sin Ψ⎠ 0
⎛cos sin ⎞ e′ = ⎝sin sin ⎠ cos where 0 denotes the x-axis nominal antenna ofset, as given in table 1. By inverting the sign convention in (5) for simplicity, the residual signal can be expressed as: 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 = [︀ ˚() ˚()]︀
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 = p · e′ = 0 (cos Ψ cos + sin Ψ sin ) sin = 0 cos Ψ · sin cos + 0 sin Ψ · sin sin 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): = [︀ 1() 2() . . . ()]︀ (9) ˚ = 0 · cos Ψ , ˚ = 0 · sin Ψ , Ψ = 2 (˚, ˚)
In this section, a set of estimation techniques will be presented, which solve diferent 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.
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: (6) (7) (8) (10) 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: ⎛ 1
˚ A = ⎜⎜⎜ .˚2
. ⎜ . ⎝ ˚ ˚1 ⎞ ⎛ sin 1 cos 1 sin 1 sin 1 ⎞ ˚2 ⎟⎟⎟ = ⎜⎜⎜ sin 2 cos 2 sin 2 sin 2 ⎟ ...˚ ⎠⎟ ⎝sin ...cos sin 1 s...in ⎟⎟⎠ 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 ·
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 [
‖Ap − ‖22
Cp = d
(11)
(12)
(13)
(14)
(15)
(16)
where Cp = d is a linear constraint. Then the system of equations introduced above can be solved
using the Karush-Kuhn-Tucker equations [
˚2 + ˚2 = 02 · (︀ cos2 Ψ + sin2 Ψ)︀ = 02 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 − ‖22 s.t. [︀ ˚ ˚ ]︀ · p+1 = 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 ofset, i.e., the radius of the antenna rotation circle, as presented in table 1.
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 [
= 0 cos Ψ · sin cos + 0 sin Ψ · sin sin + ˚ After applying this change, the grand design matrix and unknowns vector change to: ⎛ sin 1 cos 1 sin 1 sin 1 1⎞ A = ⎝⎜⎜⎜ sin 2 ...cos 2 sin 2 ...sin 2 1... ⎟⎠⎟⎟ , p+1 = ⎝⎛˚˚˚+++111⎠⎞
sin cos sin 1 sin 1 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 ofset bias is satellite-related and only one satellite is observed, there will only be one value per epoch to estimate for all stations. Finally, the constraint is reformulated as in (12), with C = (︀ ˚ ˚ 0︀) and remains a scalar, = [︀ 20 + ˚2 + ˚2 ]︀ /2. Then the optimization problem is formulated using the updated grand design matrix A and unknowns vector p+1 from (16).
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 diferent 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 cutof 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
efects 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 [
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 [
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 diferent 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 are estimated roughly between 0.058 °/ − 0.065 °/, while the rates for noon turns are roughly between 0.088 °/ − 0.1 °/. In what regards the rotation rate estimates, all results are consistent with values computed using RKPP in previous studies.
Raw data
Perc. (%)
Another feature that is not always properly modeled when estimating yaw values is the efect of the
bias applied to the observed Sun angle. It’s efect on attitude maneuvers was observed in [
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 [
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
4. In [17] it was stated that the inversion caused by the − 0.5° bias afects exclusively noon maneuvers,
meaning that midnight maneuvers do not present inverted rotation rates. At the same time, [
After analyzing the efects of linear drift on the rotation rates and efects 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.
Figure 9 showcases two sets of residuals captured by ≈ 30 stations each during two maneuvers: the left side figure shows a noon maneuver performed by PRN 26, while the right side figure shows a midnight maneuver performed by PRN 27. In both cases, the values in blue depict the originally estimated residual signals, while the values in red show the residuals after applying corrections. The applied corrections were computed based on the estimated yaw values and clock biases, by estimating the diferences between their real and estimated values. Afterwards, they were then modeled into residual biases based on (16) and subtracted from the original signals. The initial residual values have increased values during the maneuvers, and a noticeable bump in the residual signal can be observed. The corrected residuals present lower values during both maneuvers.
Finally, the efect 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 ifgure, 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 [
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 efects on rate inversion, as well as efects 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, afected 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.
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 rr, onto the plane normal to the nadir direction, rr⊥. If a circular orbit is assumed, as is approximately the case for GPS satellites, the following equations hold: h = r × v , ℎ = · , rr⊥ = cos ·
+ sin · ℎ v
h rr⊥ = rr − ︁( rr · r )︁ 2 · r 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 rr⊥ and v, respectively h, can be expressed: This then leads to the formula for computing the azimuth angle, , and the nadir angle, : rr⊥ · v = ⊥ · cos rr⊥ · h = ⊥ · ℎ sin = ⊥ · sin = 2
︂( rr⊥ · h , rr⊥ · v , cos = − ︂) (rr − r) · r |rr − r| · (17) (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28)
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:
˚+1 = ˚ + Δ˚ , ˚+1 = ˚ + Δ˚ 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) with the approximations from (22): with 0 ≈ 39.4 cm denoting the radius of the antenna rotation circle, as presented in table 1. 02 = ˚2+1 + ˚+1 = (˚ + Δ˚)2 + (︀ ˚ + Δ˚)︀ 2
2 02 = (︀ ˚2 + 2˚ Δ˚ + Δ˚2)︀ + (︀ ˚2 + 2˚ Δ˚ + Δ˚2)︀ Considering that the terms Δ˚2 and Δ˚2 become very close to zero, the constraint is simplified to: 2˚ Δ˚ + 2˚ Δ˚ = 02 − (︀ ˚2 + ˚2 )︀ And finally, by replacing rewritten as:
Δ˚ = ˚+1 − ˚ and respectively Δ˚ = ˚+1 − ˚ , the constraint can be 2˚ (˚+1 − ˚ ) + 2˚
︀( ˚+1 − ˚ )︀ = 02 − (︀ ˚2 + ˚2 )︀ 2˚ ˚+1 + 2˚
˚+1 = 02 + (︀ ˚2 + ˚2 )︀ ︀( ˚ ˚ )︀ · ︂( ˚+1︂) ˚+1 = 1 2 where ˚ and ˚ denote the known variables, estimated for the previous iteration. The final form of the constraint will then become: [17] J. Kouba, Noon turns for deep eclipsing block iif gps satellites, note prepared for IGS ACs (2013). [18] R. Dach, S. Schaer, U. Hugentobler, M. Meindl, G. Beutler, Gnss satellite clock estimation at code, IGS Workshop 2006, Darmstadt (2006). URL: https://www.bernese.unibe.ch/publist/2006/post/rd_ igsws_06_05.pdf.