=Paper=
{{Paper
|id=Vol-3719/107
|storemode=property
|title=Estimation of GPS Block IIF Attitude Using Residual Signals
|pdfUrl=https://ceur-ws.org/Vol-3719/paper2.pdf
|volume=Vol-3719
|authors=Alexandra Reitu,Urs Hugentobler,Bingbing Duan,Kaan Çelikbilek,Elena Simona Lohan,Thodoris Spanos,Fran Fabra,José A. López-Salcedo,Gonzalo Seco-Granados,Nikos Kanistras,Ivan Lapin,Vassilis Paliouras,Izadora A. Ramos,Vinícius M. G. B. Cavalcanti,Felipe O. Silva,Danilo A. de Lima
|dblpUrl=https://dblp.org/rec/conf/wiphal/ReituHD24
}}
==Estimation of GPS Block IIF Attitude Using Residual Signals==
https://ceur-ws.org/Vol-3719/paper2.pdf
Estimation of GPS Block IIF Attitude Using Residual
Signals
Alexandra Reitu1,* , Urs Hugentobler1 and Bingbing Duan1
1
Dept. of Aerospace and Geodesy, Technical University of Munich, Germany
Abstract
The scope of this work is to present a new approach for estimating the attitude of GNSS satellites with horizontal
offset 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.
Keywords
Reverse Kinematic Point Positioning, Attitude Estimation, Residual Signals
1. 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].
1.1. 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]:
− tan 𝛽
Ψ𝑛 = 𝐴𝑇 𝐴𝑁 2 (1)
sin 𝜇
WIPHAL 2024: Work-in-Progress in Hardware and Software for Location Computation, June 25-27, 2024, Antwerp, Belgium
*
Corresponding author.
$ alexandra.reitu@tum.de (A. Reitu); urs.hugentobler@tum.de (U. Hugentobler); bingbing.duan@tum.de (B. Duan)
https://github.com/reitualexandra (A. Reitu)
0009-0005-2498-2381 (A. Reitu); 0000-0003-0801-8259 (U. Hugentobler); 0000-0002-0143-835X (B. Duan)
© 2024 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR
ceur-ws.org
Workshop ISSN 1613-0073
Proceedings
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 spacecraft-
fixed 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).
y x
z
µ ß
Midnight Noon
Figure 1: Example of satllite orbit depicting the noon and midnight points, alongside the 𝛽 and 𝜇 angles
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]:
(︂ )︂
− tan 𝛽
Ψ𝑎 = 𝐴𝑇 𝐴𝑁 2 + 𝐵(𝑏, 𝛽, 𝜇) (2)
sin 𝜇
(︂ )︂
0.0175 · 𝑏
𝐵(𝑏, 𝛽, 𝜇) = 𝐴𝑆𝐼𝑁 , cos 𝐸 = cos 𝛽 · sin 𝜇 (3)
sin 𝐸
where 𝑏 = −0.5°, and 𝐵 denotes the actual yaw bias introduced by 𝑏.
1.2. 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 ionosphere-
free 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°
bias 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.
2. 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|.
2.1. 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 − 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 offset 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 |
Figure 2: Yaw angle and it’s dependence on azimuth and nadir
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 − rr
𝜌 = |r + p − rr | − |r − rr | = · p = e · p = −e′ · p (5)
|r − rr |
where e = |r−r
r−rr
r|
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:
⎛ ⎞ ⎛ ⎞
cos Ψ cos 𝛼 sin 𝜂
p = 𝑥0 ⎝ sin Ψ ⎠ e′ = ⎝ sin 𝛼 sin 𝜂 ⎠ (6)
0 cos 𝜂
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)
3. 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.
3.1. 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:
]︀𝑇
(10)
[︀
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
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 𝜕𝜌1
⎞ ⎛ ⎞
𝜕𝑥
˚ 𝜕𝑦
˚ sin 𝜂1 cos 𝛼1 sin 𝜂1 sin 𝛼1
⎜ 𝜕𝜌2 𝜕𝜌2 ⎟ ⎜ sin 𝜂2 cos 𝛼2
⎜ 𝜕𝑥
˚ 𝜕𝑦
˚ ⎟ sin 𝜂2 sin 𝛼2 ⎟
⎜ ..
A=⎜ .. ⎟
⎟=⎜ .. .. (11)
⎜ ⎟
. .
⎟
⎝ . . ⎠ ⎝ ⎠
𝜕𝜌𝐾 𝜕𝜌𝐾 sin 𝜂𝐾 cos 𝛼𝐾 sin 𝜂1 sin 𝛼𝐾
𝜕𝑥
˚ 𝜕𝑦
˚
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
(︀ 𝑇 )︀−1 𝑇
pseudo-inverse: p = A A A ·𝜌
3.2. Constraining antenna offset
The previous section described a simple model for estimating the yaw angle using only acquired obser-
vations 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 − 𝜌‖22
𝑝
(12)
s.t. Cp = d
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 = 𝑥20 · cos2 Ψ + sin2 Ψ = 𝑥20 (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 − 𝜌‖22
𝑝
1 (︀ 2 (14)
s.t. 𝑥2𝑗 + ˚
𝑦 2𝑗
[︀ ]︀ )︀
˚ 𝑦 𝑗 · p𝑗+1 =
𝑥𝑗 ˚ 𝑥0 + ˚
2
]︀𝑇
where p𝑗+1 = ˚ 𝑦 𝑗+1 denotes the solution computed for time epoch 𝑡 at iteration 𝑗 + 1. The
[︀
𝑥𝑗+1 ˚
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.
3.3. Estimation of satellite clock errors
Furthermore, it can be taken into account that some of the errors caused by poorly modeled attitude
˚𝑠 , and are no longer included in the
were absorbed by satellite clock corrections into the bias term 𝛿𝜏
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:
⎛ ⎞
sin 𝜂1 cos 𝛼1 sin 𝜂1 sin 𝛼1 1 ⎛ ⎞
⎜ sin 𝜂2 cos 𝛼2 sin 𝜂2 sin 𝛼2 1⎟ ˚
𝑥𝑗+1
A=⎜ .. .. .. ⎟ , p𝑗+1 = ⎝ ˚𝑦 𝑗+1 ⎠ (16)
⎜ ⎟
⎝ . . .⎠ ˚ 𝑠
𝑐𝛿 𝜏𝑗+1
sin 𝜂𝐾 cos 𝛼𝐾 sin 𝜂1 sin 𝛼𝐾 1
Figure 3: Examples of raw residual signals for two noon turns of PRN 26
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.
Finally, the constraint is reformulated as in (12), with C = ˚ 𝑥𝑗 ˚𝑦 𝑗 0 and 𝑑 remains a scalar,
𝑑 = 𝑥20 + ˚ 𝑦 2𝑗 /2. Then the optimization problem is formulated using the updated grand design
𝑥2𝑗 + ˚
[︀ ]︀
matrix A and unknowns vector p𝑗+1 from (16).
4. 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].
4.1. 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.
Figure 4: Estimated yaw angle for a set of noon turns, performed by PRN 26 between days 352 - 356 of year 2020
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
Figure 5: Estimated yaw angle for a set of midnight turns, performed by PRN 25, 27 and 32
Figure 6: Estimated (blue) and nominal yaw rates (red) as function of 𝛽, for filtered data
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.
Noon rate (°/𝑠) Midn rate (°/𝑠)
Perc. (%)
Avg. Med. Avg. Med.
100 0.0864 0.0863 0.0665 0.0647
Raw
95 0.0885 0.0872 0.0659 0.0647
data 67 0.0848 0.0863 0.0681 0.0656
100 0.1044 0.1044 0.0579 0.0579
Filtered
95 0.1039 0.1044 0.0572 0.0579
data 67 0.1042 0.1044 0.0583 0.0579
Table 2
Average and median yaw (rotation) rate values for different data sets and for different percentiles
Figure 7: Estimated yaw rates for noon (red) and midnight turns (blue), as function of 𝛽, for filtered data
4.2. 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.
Figure 8: Estimated yaw angle for a set of inverted turns performed by PRN 25 (noon) and 26 (midnight) on
days 355 and 358 of 2020
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
Maneuver PRN 𝛽𝑚𝑖𝑛 (°) 𝛽𝑚𝑎𝑥 (°) 𝛽 (°) Sign (±)
N1 27 -0.01 -0.08 -0.045 -
N1 25 -0.03 -0.13 -0.08 -
N2 25 -0.44 -0.54 -0.49 -
N1 26 -0.61 -0.71 -0.66 +
N1 25 -0.85 -0.95 -0.9 -
N2 27 -0.88 -0.95 -0.915 +
Table 3
Noon maneuvers observed for −1° ≤ 𝛽 ≤ 0°
Maneuver PRN 𝛽𝑚𝑖𝑛 (°) 𝛽𝑚𝑎𝑥 (°) 𝛽 (°) Sign (±)
M1 26 -0.1 0 -0.05 -
M1 27 -0.15 -0.22 -0.185 +
M1 25 -0.23 -0.33 -0.28 +
M2 27 -0.45 -0.52 -0.485 +
M1 26 -0.81 -0.91 -0.86 +
Table 4
Midnight maneuvers observed for −1° ≤ 𝛽 ≤ 0°
4. 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°.
4.3. 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.
Figure 9: Examples of original (blue) and corrected residuals (red) for two turn
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 differences 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.
4.4. 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].
Figure 10: Detrended clock products before and after correcting them using estimated clock biases, PRN 25,
days 354 and 355 of year 2020
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.
5. 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.
6. 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 instanta-
neous 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:
v h
h = r × v , ℎ = 𝑟 · 𝑣 , rr⊥ = 𝑟𝑟 cos 𝛼 · + 𝑟𝑟 sin 𝛼 · (17)
𝑣 ℎ
(︁ r · r )︁
r
rr⊥ = rr − ·r (18)
𝑟2
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:
rr⊥ · v = 𝑟𝑟 ⊥ · 𝑣 cos 𝛼 (19)
rr⊥ · h = 𝑟𝑟 ⊥ · ℎ sin 𝛼 = 𝑟𝑟 ⊥ · 𝑟𝑣 sin 𝛼 (20)
This then leads to the formula for computing the azimuth angle, 𝛼, and the nadir angle, 𝜂:
(︂ )︂
rr⊥ · h (rr − r) · r
𝛼 = 𝐴𝑇 𝐴𝑁 2 , rr⊥ · v , cos 𝜂 = − (21)
𝑟 |rr − r| · 𝑟
7. 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:
˚
𝑥𝑗+1 = ˚
𝑥𝑗 + Δ𝑥
˚, ˚
𝑦 𝑗+1 = ˚
𝑦 𝑗 + Δ𝑦
˚ (22)
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):
)︀2
𝑥20 = ˚𝑥2𝑗+1 + ˚𝑦 2𝑗+1 = (𝑥 ˚)2 + ˚ (23)
(︀
˚𝑗 + Δ𝑥 𝑦 𝑗 + Δ𝑦
˚
𝑥20 = ˚ (24)
(︀ 2
˚2 + ˚
)︀ (︀ 2
˚2
)︀
𝑥𝑗 + 2𝑥˚𝑗 Δ𝑥 ˚ + Δ𝑥 𝑦 𝑗 + 2𝑦
˚𝑗 Δ𝑦 ˚ + Δ𝑦
Considering that the terms Δ𝑥
˚2 and Δ𝑦
˚2 become very close to zero, the constraint is simplified to:
˚ = 𝑥20 − ˚ (25)
(︀ 2
𝑦 2𝑗
)︀
2𝑥
˚𝑗 Δ𝑥
˚ + 2𝑦
˚𝑗 Δ𝑦 𝑥𝑗 + ˚
And finally, by replacing Δ𝑥˚=˚ 𝑥𝑗+1 − ˚𝑥𝑗 and respectively Δ𝑦 ˚=˚ 𝑦 𝑗+1 − ˚ 𝑦 𝑗 , the constraint can be
rewritten as:
𝑦 𝑗 = 𝑥20 − ˚ (26)
(︀ 2
𝑦 2𝑗
(︀ )︀ )︀
2𝑥 ˚𝑗+1 − ˚
˚𝑗 (𝑥 𝑥𝑗 ) + 2𝑦
˚𝑗 ˚ 𝑦 𝑗+1 − ˚ 𝑥𝑗 + ˚
𝑦 𝑗+1 = 𝑥20 + ˚ 𝑥2𝑗 + ˚
𝑦 2𝑗 (27)
(︀ )︀
2𝑥
˚𝑗˚𝑥𝑗+1 + 2𝑦˚𝑗˚
where ˚
𝑥𝑗 and ˚𝑦 𝑗 denote the known variables, estimated for the previous iteration. The final form of the
constraint will then become:
(︂ )︂
)︀ ˚ 𝑥𝑗+1 1 (︀ 2
𝑥2𝑗 + ˚
𝑦 2𝑗 (28)
(︀ )︀
˚ 𝑦𝑗 ·
𝑥𝑗 ˚ = 𝑥0 + ˚
˚𝑦 𝑗+1 2
with 𝑥0 ≈ 39.4 cm denoting the radius of the antenna rotation circle, as presented in table 1.
References
[1] X. Cao, S. Zhang, K. Kuang, T. Liu, K. Gao, The impact of eclipsing gnss satellites on the precise
point positioning, Remote Sensing 10 (2018) 94. URL: https://www.researchgate.net/publication/
322408094_The_Impact_of_Eclipsing_GNSS_Satellites_on_the_Precise_Point_Positioning. doi:10.
3390/rs10010094.
[2] S. Strasser, S. Banville, A. Kvas, S. Loyer, T. Mayer-Gürr, Comparison and generalization of gnss
satellite attitude models, EGU General Assembly (2021). URL: https://graz.pure.elsevier.com/
en/publications/comparison-and-generalization-of-gnss-satellite-attitude-models. doi:10.1002/
andp.19053221004.
[3] O. Montenbruck, R. Schmid, F. Mercier, P. Steigenberger, C. Noll, R. Fatkulin, S. Kogure, A. S.
Ganeshan, Gnss satellite geometry and attitude models, Advances in Space Research 56 (2015)
1015–1029. URL: https://www.sciencedirect.com/science/article/pii/S0273117715004378. doi:10.
1016/j.asr.2015.06.019.
[4] C. Rodriguez-Solano, U. Hugentobler, P. Steigenberger, S. Lutz, Impact of earth radiation pressure
on gps position estimates, Journal of Geodesy 86 (2012) 309–317. URL: https://link.springer.com/
article/10.1007/s00190-011-0517-4. doi:10.1007/s00190-011-0517-4.
[5] Y. E. Bar-Sever, A new model for gps yaw attitude, Journal of Geodesy 70 (1996) 714–723. URL:
https://link.springer.com/article/10.1007/BF00867149. doi:10.1007/BF00867149.
[6] D. Kuang, S. Desai, A. Sibois, Observed features of gps block iif satellite yaw maneuvers
and corresponding modeling, GPS Solutions (2017). URL: https://www.researchgate.net/
publication/306924379_Observed_features_of_GPS_Block_IIF_satellite_yaw_maneuvers_and_
corresponding_modeling. doi:0.1007/s10291-016-0562-9.
[7] F. Dilssner, T. Springer, W. Enderle, Gps iif yaw attitude control during eclipse season, AGU Fall
Meeting, San Francisco (2011). URL: http://acc.igs.org/orbits/yaw-IIF_ESOC_agu11.pdf.
[8] F. Dilssner, Gps iif-1 satellite antenna phase center and attitude modeling, InsideGNSS (2010).
URL: https://www.insidegnss.com/auto/sep10-Dilssner.pdf.
[9] A. Sibois, A. Sibthorpe, A multi-year reanalysis of gps block ii/iia and iif satellite yaw maneuvers
by means of reverse kinematic point positioning technique, IGS Workshop 2018, Wuhan
(2018). URL: https://www.researchgate.net/publication/328942357_A_multi-year_reanalysis_of_
GPS_Block_IIIIA_and_IIF_satellite_yaw_maneuvers_by_means_of_reverse_kinematic_point_
positioning_technique.
[10] O. Colombo, T. Thomas, D. D. Rowlands, S. B. Luthcke, Testing a reverse kinematic point positioning
technique for possible operational use in gps data analysis at nasa’s goddard space flight center,
IGS Workshop 2016, Sydney (2016). URL: https://science.gsfc.nasa.gov/earth/geodesy/content/
uploadFiles/highlight_files/IGS_Sydney.Poster_OscarColombo_etal_.pdf.
[11] R. Dach, S. Lutz, P. Walser, P. Fridez, Bernese gnss software version 5.2., Astronomical Institute,
University of Bern, Bern Open Publishing (2015). URL: http://ftp.aiub.unibe.ch/BERN52/DOCU/
DOCU52.pdf. doi:10.7892/boris.72297.
[12] S. Lutz, G. Beutler, S. Schaer, Code’s new ultra-rapid orbit and erp products for the igs, GPS
Solutions 20 (2016) 239–250. URL: https://link.springer.com/article/10.1007/s10291-014-0432-2.
doi:10.1007/s10291-014-0432-2.
[13] R. Dach, S. Schär, D. Arnold, E. Brockmann, M. S. Kalarus, L. Prange, P. Stebler, A. Jäggi, Code
final product series for the igs, Astronomical Institute, University of Bern (2020). doi:10.48350/
185646.
[14] R. L. Burden, J. D. Faires, Numerical Analysis, 9 ed., Brooks/Cole, 2011.
[15] J. Scott, M. Tůma, Solving large linear least squares problems with linear equality constraints,
BIT Numerical Mathematics (2022) 1572–9125. URL: https://link.springer.com/article/10.1007/
s10543-022-00930-2. doi:10.1007/s10543-022-00930-2.
[16] D. Kuang, S. Desai, A. Sibois, Observed features of gps block iif satellite yaw maneuvers and
corresponding modeling, GPS Solutions 21 (2017) 739–745. URL: https://link.springer.com/article/
10.1007/s10291-016-0562-9. doi:10.1007/s10291-016-0562-9.
[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.