# Performance Trades for Multiantenna GNSS Multisensor Attitude Determination Systems

Authors: Nathan A. Tehrani; Jason N. Gross
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

We present various performance trades for multiantenna global navigation satellite system (GNSS) multisensor attitude estimation systems. In particular, attitude estimation performance sensitivity to various error sources and system configurations is assessed. This study is motivated by the need for system designers, scientists, and engineers of airborne astronomical and remote sensing platforms to better determine which system configuration is most suitable for their specific application. In order to assess performance trade-offs, the attitude estimation performance of various approaches is tested using a simulation that is based on a stratospheric balloon platform. For GNSS errors, attention is focused on multipath, receiver measurement noise, and carrier-phase breaks. For the remaining attitude sensors, different performance grades of sensors are assessed. Through a Monte Carlo simulation, it is shown that, under typical conditions, sub-0.1-degree attitude accuracy is available when using multiple antenna GNSS data only, but that this accuracy can degrade to degree level in some environments warranting the inclusion of additional attitude sensors to maintain the desired level of accuracy. Further, we show that integrating inertial sensors is more valuable whenever accurate pitch and roll estimates are critical.


---

## Body

## 1. Introduction

For any airborne sensing platform, the pointing accuracy is dependent on and can be limited to the accuracy of the onboard attitude solution [1, 2]. As such, a key to high pointing accuracy is a robust attitude-determination system. This paper outlines the development, simulation, and testing of a multisensor attitude determination algorithm intended for airborne astronomical and remote sensing platforms.Attitude determination using multiantenna global navigation satellite system (GNSS) observations is a well-established technology, first proposed by Cohen and Parkinson in 1991 for spacecraft applications [3]. It was also adapted for aircraft use [4] and tested by the same author [5]. Multiantenna GNSS attitude determination has been tested on ground, waterborne, and flight vehicles [6], and the technology has matured to the level that multiple commercially available products [7, 8] are available. The technology has been used for remote sensing platforms since shortly after its proposal [9], and it is in use on multiple stratospheric balloon platforms [10].There has been considerable effort to simulate gyroscope-free attitude determination using 3-axis magnetometers, 2-axis sun sensors, or both, for spacecraft applications [11]. Highlights include the use of a magnetometer-only sun-pointing algorithm by Ahn and Lee in 2003 [12]. The magnetometer-derived attitude was within 3° of a gyroscope-derived “truth.” Psiaki modeled an orbit- and attitude-determination algorithm [13]. Using a 10 nT 3-axis magnetometer and a 0.005°-σ sun sensor, this approach showed less than 0.1° error in all axes. Crassidis and Markley created a sun sensor and magnetometer Kalman filter and showed that a magnetometer-only attitude estimate is markedly improved (error reduced by approximately half) with the inclusion of sun sensor data [11]. Considering the increasing presence of small uninhabited aerial system (UAS) applications in remote sensing, sensor fusion promises high attitude performance with low-cost, lightweight, and compact hardware [10, 14].This paper outlines the design of a GNSS-based attitude estimator that is then optionally augmented with various other attitude sensors for use in other airborne remote sensing or astronomy applications and offers the contribution of assessing performance sensitivity to various design configuration and common error sources. This paper makes up part of the first author’s graduate thesis [15] and is a significant extension upon our previous conference paper [16] in which we improve upon the presented algorithm formulation and revamp our Monte Carlo simulation study design. Through a simulation that is built upon multiple sensor models, the GNSS-only-based attitude solution is shown to work well, but is significantly improved when additional sensors are optimally fused. Specifically, we show the benefits of including inertial, sun sensor, and magnetometer measurements. This paper is expected to aid system designers and scientific investigators as they propose or implement an attitude determination system that is required to meet a certain accuracy level.The remainder of this paper is organized as follows. In Section2, the design of the baseline estimation filter and attitude estimation filter is discussed. In Section 3, the simulation environment and the sensor data simulation are discussed. In Section 4, the performance of the GNSS-based and multisensor attitude estimators is presented and discussed. Section 5 summarizes this study’s findings and discusses future work.


## 2. Attitude Estimation

### 2.1. Algorithm Overview

Figure1 shows the overall algorithm considered in this study. First, a carrier-phase differential GNSS filter, as detailed in Section 2.2, estimates the baseline separations between antennas. Next, this information is used within a GNSS-only multiple antenna attitude estimator as described in Section 2.3. Finally, the resulting estimated attitude state is optionally fused with a multisensor estimator that also incorporates information from inertial sensors, magnetometers, and sun sensors, as discussed in Section 2.4.Figure 1
Block diagram showing the three main estimators: antenna baseline filter, GNSS-only attitude estimator, and multisensor attitude estimator. This figure is an updated version of what was presented in the conference paper version of this work [16].


### 2.2. Antenna Baseline Filter

This attitude estimation algorithm begins with antenna baseline determination using a Kalman filter. The Kalman filter estimator is a linear state estimator which was developed by Kalman, Swerling, and Bucy in a series of papers which detail its formulation and implementation [17–19]. A full derivation of the Kalman filter can be found in Crassidis and Markley and Groves [11, 20]. This section will focus on the state, covariance, and tuning of the Kalman filter applied in this study.This implementation uses pseudorange and carrier-phase differential GNSS (CD-GNSS) measurements to estimate the relative position vectors between the antennas [21]. The state vector, x, for this filter consists of the relative position vector components between antennas A and B, xA,B, yA,B, and zA,B, and a set of double-differenced pseudoranges and carrier-phase biases, NA,B.
(1)xbaseline=xA,ByA,BzA,BNA,B1,k⋮NA,Bj,k.The measurement models used to model the double-differenced carrier-phase observables follow the same approach outlined in [22], as is discussed next.The model for an undifferenced GNSS pseudorange,ρ, and carrier-phase measurement, ϕ, (with units of carrier cycles) is given as [21]
(2)ρ=r+Iρ+Tρ+cδtu−δts+ερ,ϕ=λ−1r+Iϕ+Tϕ+cλδtu−δts+N+εϕ,where λ is the wavelength corresponding to the frequencies L1 and L2 and expressed in meters. The geometric range r between the receiver and GNSS satellite is also expressed in meters, as are the ionospheric and tropospheric delays I and T. The speed of light c is expressed in meters per second. The clock biases of the receiver and satellite, δtu and δts, respectively, are expressed in seconds. The unmodeled error sources, which include multipath and thermal noise, are included in ε in units of meters.First, range and phase measurements for the master antennaA (antenna 1) and B (antennas 2, 3, or 4) are differenced to form single-differenced phase measurements:
(3)ΔρA,Bj=rA,Bj+cδtA,B+ερ,A,Bj,ΔϕA,Bj=λ−1rA,Bj+cλδtA,B+NA,Bj+εϕ,A,Bj.Within (3), due to the very short baseline separation between the antennas, the atmospheric delays completely cancel along with any satellite clock bias and ephemeris errors. Next, the single differenced measurements are then differenced between satellites. For example, between satellite j and reference satellite k,
(4)∇Δρf,A,Bj,k=rA/Bk∣k−1+ερ,A,Bj,k,∇Δϕf,A,Bj,k=−λ−11¯Aj−1¯AkTrA/Bk∣k−1+NA,Bj,k+εϕ,A,Bj,k,where the remaining receiver clock bias errors are eliminated, leaving only the unknown phase biases NA,Bj,k, which are known to be integers.Because the GPS and GLONASS satellite constellations operate at different frequencies, both a GPS and a GLONASS satellite are used as separate reference satellites [23]. GLONASS satellites operate using frequency division multiple access (FDMA), and the wavelength varies from satellite to satellite. The resulting interchannel bias is negligible when using like receivers, as in this model [24]. Because of this, double differences were only formed within each satellite constellation (i.e., GPS and GLONASS) and a different reference satellite was identified for each.The observation matrix,H, transforms the state x to predicted measurements. The first three rows of this filter’s observation matrix consist of the difference of two three-component unit vectors, which point from the user to the satellite used in the double differencing and to a reference satellite, R.
(5)H=ux1−uxRuy1−uyRuz1−uzRλL1,1⋮⋮⋮⋱uxn−uxRuyn−uyRuzn−uzRλL1,nux1−uxRuy1−uyRuz1−uzRλL2,1⋮⋮⋮⋱uxn−uxRuyn−uyRuzn−uzRλL2,n.The measurement vector,z, consists of double-differenced pseudorange and phase measurements for each satellite relative to the reference satellite, including measurements for each L1 and L2 frequencies:
(6)z=∇ΔρL1,A,Bi,k∇ΔρL2,A,Bi,k∇ΔϕL1,A,Bi,k∇ΔϕL2,A,Bi,k.Operating in parallel with the baseline estimation Kalman filter, the floating point estimated phase biases (for GPS satellites only);NA,Bj,k and their estimated error-covariance are fed into an integer ambiguity resolution algorithm. In particular, the least-squares ambiguity decorrelation adjustment (LAMBDA) method [25] is used to determine the integer biases and adjust the estimated relative positions.The Kalman filter process noise,Q, and measurement noise, R, and initial error-covariance, P0, assumptions selected for the differential GNSS baseline estimator are outlined in Table 1:Table 1
Baseline filter assumed parameters, same as in [16].
Filter parameter
Assumed values
State error covarianceP0
Baseline states: 1 mAmbiguity states: 225 m
Measurement noise covarianceR
σ
ϕ
=
4
⋅
10
−
4
m
Process noise covarianceQ
Attitude states:in‐runbias⋅10−2m/sAmbiguitystates:0m/s


### 2.3. Baseline to Attitude

Once the antenna relative baselines with respect to a master antenna are estimated using the baseline estimation filter, an earth-centered, earth-fixed (ECEF) antenna relative position matrix,Le, is constructed at each epoch by vertically concatenating the estimate relative vectors of each of nonmaster antenna, as adopted from Cohen and Parkinson [3]:
(7)Le=x2,ey2,ez2,ex3,ey3,ez3,ex4,ey4,ez4,e.This matrix is used to estimate the platform attitude given by the antenna baseline vectors. The known bodycentric antenna coordinate matrix,Lb, with its origin defined as the reference antenna’s position, makes up the list of reference vectors.Next, the singular value decomposition (SVD) method, as proposed by Markley and Mortari and described in (9), (10), and (11), is used to find the rotation matrix Ceb between the ECEF, e, and body, b, frames [26]. This is then transformed into the local navigation frame using, n, the earth-to-nav rotation, which is dependent on the master antennas’ latitude and longitude and is defined in many texts [20]:
(8)Cnb=CneCeb.This solution requires the construction of a matrixB using the measured vectors vi and reference vectors wi:
(9)B=∑i=1nviwiT.A SVD is performed onB, resulting in unitary matrices U and V. A diagonal, 3 × 3 matrix M is constructed using the determinants of U and V:
(10)M=11UV.By multiplyingU, M, and V, one can find the rotation matrix R:
(11)R=UMVT.The SVD rotation solution also provides a straightforward attitude error covariance matrix which was used as an error covariance matrix in the GNSS-only attitude filter and for the measurement covariance of the GNSS-attitude estimates when fusing the solution in a multisensor Kalman filter [27].In order to obtain the attitude error covariance matrix, the matrixB is multiplied by the transpose of the nontransformed rotation matrix Ceb.
(12)D=BCebT.D can then be used to find the inverse of the error covariance matrix:
(13)P−1=trD∗I3×3−D.This is then scaled by multiplying the nominal standard deviation of attitude error attributed to noise in the baseline measurements, in this case 0.01°, to obtain a measurement error covariance matrix to be used for GNSS measurements in the nonlinear filter.The attitude dilution of precision, as proposed by Yoon and Lundberg, and will be used in the paper to assess trades with respect to multiple GNSS constellations, is a similar metric which assesses the ability to measure the Euler angles [28] depending on the GNSS satellite constellation geometry. It is defined as [28]
(14)ADOP=trnI−SST−1,where n is the number of satellites in view, I is the 3 × 3 identity matrix, and S is a 3 × N matrix comprising the unit vectors to each satellite, including the reference satellite [28].


### 2.4. Multisensor Unscented Kalman Filter

Finally, a nonlinear estimator is used for attitude determination using all sensor data, due to the nonlinear nature of attitude to sensor observation transformation. Several nonlinear attitude determination methods exist. The QUEST, or quaternion estimation method, seeks the unique quaternion solution for a set of vector measurements and reference vectors [29]. The RE-QUEST, or recursive QUEST, applies the same method, but rather than solving for a single epoch, it uses a filtering approach to find the time-varying attitude profile [30]. Other filtering techniques include the multiplicative extended Kalman filter (M-EKF) [31] and the quadratic extended Kalman filter (Q-EKF) [32].In this paper, an unscented Kalman filter (UKF) was chosen for its nonlinear transformation ability. The UKF has been shown to perform nearly identically an EKF for GPS/INS attitude estimation applications [33, 34], but offers the implementation benefit of not requiring the analytic evaluation of linearized partial derivatives of the process and observation equations. The details of the UKF implementation followed in this study are offered in the tutorial paper by Rhudy et al. [34], and as such, these details are not discussed in detail herein. In this paper, an outline of the state vector, state prediction fx, and observation functions hx for each measurement update is discussed.This unscented Kalman filter utilizes inertial measurements for its state prediction, with bias estimation. The GNSS Euler angles, as well as magnetometer and sun sensor data, are used as measurement updates. Magnetometer measurement updates occur at each filter time step, which take place at 50 Hz intervals. GNSS and sun sensor updates occur at 10 Hz intervals.The state vector,x estimated in the multisensor filter is given as
(15)xUKF=ϕθψbpbqbr.where ϕ, θ, and ϕ are the platform’s roll, pitch, and yaw and bp,q,r is the time-varying biases of the IMU’s roll rate, p, pitch rate, q, and yaw rate, r, gyroscopes.Within the UKF framework, at each epoch, the state vector is expanded into a group of2L+1 sigma points, χ, where L=6 is the length of the estimated state vector. For each group of sigma points l, the attitude states are predicted by integrating the IMU gyro data through the attitude kinematic equations [35]:
(16)fϕ,θ,ψ:ϕiθiϕi=ϕi−1θi−1ϕi−1+1tθi−1sϕi−1tθi−1cϕi−10cϕi−1−sϕi−1cθi−10sϕcθi−1cϕi−1cθi−1pqr−bpbqbrΔt,where s⋅ represents sine, c⋅ represents cosine, and t⋅ represents tangent.The delta anglesp, q, and r, shown in Figure 2, are the measured delta angles from the gyroscope, corrected for the craft- and earth-rate rotations which were included in the IMU model.
(17)pqr=p′q′r′−pcqcrc+peqere.Figure 2
The roll, pitch, and yaw delta anglesp, q, and r.The earth-rate rotationspe, qe, and re can be found using the earth’s rotation transformed to the body frame, then converting to Euler angles [20]
(18)Δθe=ΩE0sinL0−sinL0−cosL0cosL0Cnb,where ΩE is earth’s rotation rate and L is the craft’s latitude. The craft-rate rotation component can be found in a similar way [20]
(19)Δθc=0−ωen,znωen,ynωen,zn0−ωen,xn−ωen,ynωen,xn0,ωenn=veb,EnREL+h−veb,NnRNL+h−veb,EntanLREL+h,where veb,Nn and veb,En are the north and east components of the craft’s velocity, h is the craft’s altitude, and RE is earth’s radius.Furthermore,ϕi−1, θi−1, and ψi−1 are the previous epoch’s roll, pitch, and yaw sigma points and are the first three elements of each column of χ, and bp,q,r is the sigma points corresponding to the estimated IMU bias, which is predicted as random walk parameters.
(20)fbp,q,r:bp,ibq,ibr,i=bp,i−1bq,i−1br,i−1+wb,pwb,qwb,r.The measurement-prediction matrixΨ is populated by the predicted measurement vectors using each set of sigma-points in χ. Because measurements occur at different rates in this filter, it is necessary to have different measurement updates occur at different rates. For epochs coinciding with the sun sensor and GNSS attitude measurements, each column Ψi is as follows:
(21)Ψi=Bb,xBb,yBb,z∠X∠Yϕ′θ′ψ′,where Bb, ∠X, and ∠Y are predicted magnetometer and sun sensor measurements based on the ith sigma point.The observation equations useC^nb, the direction-cosine representation of the predicted attitude states ϕ^, θ^, and ψ^. Modeling of the magnetometer measurements uses the attitude state multiplied by the earth’s magnetic field vector:
(22)hBϕ,θ,ψ:B→b=C^nbB→n.Sun sensor measurements begin with the navigation-frame solar incidence vectorVsun,n, which are multiplied by the attitude state to find the body-frame incidence vector
(23)Vsun,b=C^nbVsun,n,and the solar incidence angles ∠X and ∠Y are then calculated as follows:
(24)h∠Xϕ,θ,ψ:∠X=π2+atan2Vsun,z,bVsun,x,b,h∠Yϕ,θ,ψ:∠Y=π2+atan2sunz,bsuny,b,where atan2 is the four-quadrant tangent inverse.The measurement update matrixz consists of the simulated sensor measurement at each filter epoch. These are similar in form to the columns of Ψ:
(25)z=Bb,xBb,yBb,z∠X∠YϕGNSSθGNSSψGNSS.For epochs without GNSS and sun sensor measurements,Ψ and z consist only of magnetometer predictions and measurements.


## 2.1. Algorithm Overview

Figure1 shows the overall algorithm considered in this study. First, a carrier-phase differential GNSS filter, as detailed in Section 2.2, estimates the baseline separations between antennas. Next, this information is used within a GNSS-only multiple antenna attitude estimator as described in Section 2.3. Finally, the resulting estimated attitude state is optionally fused with a multisensor estimator that also incorporates information from inertial sensors, magnetometers, and sun sensors, as discussed in Section 2.4.Figure 1
Block diagram showing the three main estimators: antenna baseline filter, GNSS-only attitude estimator, and multisensor attitude estimator. This figure is an updated version of what was presented in the conference paper version of this work [16].


## 2.2. Antenna Baseline Filter

This attitude estimation algorithm begins with antenna baseline determination using a Kalman filter. The Kalman filter estimator is a linear state estimator which was developed by Kalman, Swerling, and Bucy in a series of papers which detail its formulation and implementation [17–19]. A full derivation of the Kalman filter can be found in Crassidis and Markley and Groves [11, 20]. This section will focus on the state, covariance, and tuning of the Kalman filter applied in this study.This implementation uses pseudorange and carrier-phase differential GNSS (CD-GNSS) measurements to estimate the relative position vectors between the antennas [21]. The state vector, x, for this filter consists of the relative position vector components between antennas A and B, xA,B, yA,B, and zA,B, and a set of double-differenced pseudoranges and carrier-phase biases, NA,B.
(1)xbaseline=xA,ByA,BzA,BNA,B1,k⋮NA,Bj,k.The measurement models used to model the double-differenced carrier-phase observables follow the same approach outlined in [22], as is discussed next.The model for an undifferenced GNSS pseudorange,ρ, and carrier-phase measurement, ϕ, (with units of carrier cycles) is given as [21]
(2)ρ=r+Iρ+Tρ+cδtu−δts+ερ,ϕ=λ−1r+Iϕ+Tϕ+cλδtu−δts+N+εϕ,where λ is the wavelength corresponding to the frequencies L1 and L2 and expressed in meters. The geometric range r between the receiver and GNSS satellite is also expressed in meters, as are the ionospheric and tropospheric delays I and T. The speed of light c is expressed in meters per second. The clock biases of the receiver and satellite, δtu and δts, respectively, are expressed in seconds. The unmodeled error sources, which include multipath and thermal noise, are included in ε in units of meters.First, range and phase measurements for the master antennaA (antenna 1) and B (antennas 2, 3, or 4) are differenced to form single-differenced phase measurements:
(3)ΔρA,Bj=rA,Bj+cδtA,B+ερ,A,Bj,ΔϕA,Bj=λ−1rA,Bj+cλδtA,B+NA,Bj+εϕ,A,Bj.Within (3), due to the very short baseline separation between the antennas, the atmospheric delays completely cancel along with any satellite clock bias and ephemeris errors. Next, the single differenced measurements are then differenced between satellites. For example, between satellite j and reference satellite k,
(4)∇Δρf,A,Bj,k=rA/Bk∣k−1+ερ,A,Bj,k,∇Δϕf,A,Bj,k=−λ−11¯Aj−1¯AkTrA/Bk∣k−1+NA,Bj,k+εϕ,A,Bj,k,where the remaining receiver clock bias errors are eliminated, leaving only the unknown phase biases NA,Bj,k, which are known to be integers.Because the GPS and GLONASS satellite constellations operate at different frequencies, both a GPS and a GLONASS satellite are used as separate reference satellites [23]. GLONASS satellites operate using frequency division multiple access (FDMA), and the wavelength varies from satellite to satellite. The resulting interchannel bias is negligible when using like receivers, as in this model [24]. Because of this, double differences were only formed within each satellite constellation (i.e., GPS and GLONASS) and a different reference satellite was identified for each.The observation matrix,H, transforms the state x to predicted measurements. The first three rows of this filter’s observation matrix consist of the difference of two three-component unit vectors, which point from the user to the satellite used in the double differencing and to a reference satellite, R.
(5)H=ux1−uxRuy1−uyRuz1−uzRλL1,1⋮⋮⋮⋱uxn−uxRuyn−uyRuzn−uzRλL1,nux1−uxRuy1−uyRuz1−uzRλL2,1⋮⋮⋮⋱uxn−uxRuyn−uyRuzn−uzRλL2,n.The measurement vector,z, consists of double-differenced pseudorange and phase measurements for each satellite relative to the reference satellite, including measurements for each L1 and L2 frequencies:
(6)z=∇ΔρL1,A,Bi,k∇ΔρL2,A,Bi,k∇ΔϕL1,A,Bi,k∇ΔϕL2,A,Bi,k.Operating in parallel with the baseline estimation Kalman filter, the floating point estimated phase biases (for GPS satellites only);NA,Bj,k and their estimated error-covariance are fed into an integer ambiguity resolution algorithm. In particular, the least-squares ambiguity decorrelation adjustment (LAMBDA) method [25] is used to determine the integer biases and adjust the estimated relative positions.The Kalman filter process noise,Q, and measurement noise, R, and initial error-covariance, P0, assumptions selected for the differential GNSS baseline estimator are outlined in Table 1:Table 1
Baseline filter assumed parameters, same as in [16].
Filter parameter
Assumed values
State error covarianceP0
Baseline states: 1 mAmbiguity states: 225 m
Measurement noise covarianceR
σ
ϕ
=
4
⋅
10
−
4
m
Process noise covarianceQ
Attitude states:in‐runbias⋅10−2m/sAmbiguitystates:0m/s


## 2.3. Baseline to Attitude

Once the antenna relative baselines with respect to a master antenna are estimated using the baseline estimation filter, an earth-centered, earth-fixed (ECEF) antenna relative position matrix,Le, is constructed at each epoch by vertically concatenating the estimate relative vectors of each of nonmaster antenna, as adopted from Cohen and Parkinson [3]:
(7)Le=x2,ey2,ez2,ex3,ey3,ez3,ex4,ey4,ez4,e.This matrix is used to estimate the platform attitude given by the antenna baseline vectors. The known bodycentric antenna coordinate matrix,Lb, with its origin defined as the reference antenna’s position, makes up the list of reference vectors.Next, the singular value decomposition (SVD) method, as proposed by Markley and Mortari and described in (9), (10), and (11), is used to find the rotation matrix Ceb between the ECEF, e, and body, b, frames [26]. This is then transformed into the local navigation frame using, n, the earth-to-nav rotation, which is dependent on the master antennas’ latitude and longitude and is defined in many texts [20]:
(8)Cnb=CneCeb.This solution requires the construction of a matrixB using the measured vectors vi and reference vectors wi:
(9)B=∑i=1nviwiT.A SVD is performed onB, resulting in unitary matrices U and V. A diagonal, 3 × 3 matrix M is constructed using the determinants of U and V:
(10)M=11UV.By multiplyingU, M, and V, one can find the rotation matrix R:
(11)R=UMVT.The SVD rotation solution also provides a straightforward attitude error covariance matrix which was used as an error covariance matrix in the GNSS-only attitude filter and for the measurement covariance of the GNSS-attitude estimates when fusing the solution in a multisensor Kalman filter [27].In order to obtain the attitude error covariance matrix, the matrixB is multiplied by the transpose of the nontransformed rotation matrix Ceb.
(12)D=BCebT.D can then be used to find the inverse of the error covariance matrix:
(13)P−1=trD∗I3×3−D.This is then scaled by multiplying the nominal standard deviation of attitude error attributed to noise in the baseline measurements, in this case 0.01°, to obtain a measurement error covariance matrix to be used for GNSS measurements in the nonlinear filter.The attitude dilution of precision, as proposed by Yoon and Lundberg, and will be used in the paper to assess trades with respect to multiple GNSS constellations, is a similar metric which assesses the ability to measure the Euler angles [28] depending on the GNSS satellite constellation geometry. It is defined as [28]
(14)ADOP=trnI−SST−1,where n is the number of satellites in view, I is the 3 × 3 identity matrix, and S is a 3 × N matrix comprising the unit vectors to each satellite, including the reference satellite [28].


## 2.4. Multisensor Unscented Kalman Filter

Finally, a nonlinear estimator is used for attitude determination using all sensor data, due to the nonlinear nature of attitude to sensor observation transformation. Several nonlinear attitude determination methods exist. The QUEST, or quaternion estimation method, seeks the unique quaternion solution for a set of vector measurements and reference vectors [29]. The RE-QUEST, or recursive QUEST, applies the same method, but rather than solving for a single epoch, it uses a filtering approach to find the time-varying attitude profile [30]. Other filtering techniques include the multiplicative extended Kalman filter (M-EKF) [31] and the quadratic extended Kalman filter (Q-EKF) [32].In this paper, an unscented Kalman filter (UKF) was chosen for its nonlinear transformation ability. The UKF has been shown to perform nearly identically an EKF for GPS/INS attitude estimation applications [33, 34], but offers the implementation benefit of not requiring the analytic evaluation of linearized partial derivatives of the process and observation equations. The details of the UKF implementation followed in this study are offered in the tutorial paper by Rhudy et al. [34], and as such, these details are not discussed in detail herein. In this paper, an outline of the state vector, state prediction fx, and observation functions hx for each measurement update is discussed.This unscented Kalman filter utilizes inertial measurements for its state prediction, with bias estimation. The GNSS Euler angles, as well as magnetometer and sun sensor data, are used as measurement updates. Magnetometer measurement updates occur at each filter time step, which take place at 50 Hz intervals. GNSS and sun sensor updates occur at 10 Hz intervals.The state vector,x estimated in the multisensor filter is given as
(15)xUKF=ϕθψbpbqbr.where ϕ, θ, and ϕ are the platform’s roll, pitch, and yaw and bp,q,r is the time-varying biases of the IMU’s roll rate, p, pitch rate, q, and yaw rate, r, gyroscopes.Within the UKF framework, at each epoch, the state vector is expanded into a group of2L+1 sigma points, χ, where L=6 is the length of the estimated state vector. For each group of sigma points l, the attitude states are predicted by integrating the IMU gyro data through the attitude kinematic equations [35]:
(16)fϕ,θ,ψ:ϕiθiϕi=ϕi−1θi−1ϕi−1+1tθi−1sϕi−1tθi−1cϕi−10cϕi−1−sϕi−1cθi−10sϕcθi−1cϕi−1cθi−1pqr−bpbqbrΔt,where s⋅ represents sine, c⋅ represents cosine, and t⋅ represents tangent.The delta anglesp, q, and r, shown in Figure 2, are the measured delta angles from the gyroscope, corrected for the craft- and earth-rate rotations which were included in the IMU model.
(17)pqr=p′q′r′−pcqcrc+peqere.Figure 2
The roll, pitch, and yaw delta anglesp, q, and r.The earth-rate rotationspe, qe, and re can be found using the earth’s rotation transformed to the body frame, then converting to Euler angles [20]
(18)Δθe=ΩE0sinL0−sinL0−cosL0cosL0Cnb,where ΩE is earth’s rotation rate and L is the craft’s latitude. The craft-rate rotation component can be found in a similar way [20]
(19)Δθc=0−ωen,znωen,ynωen,zn0−ωen,xn−ωen,ynωen,xn0,ωenn=veb,EnREL+h−veb,NnRNL+h−veb,EntanLREL+h,where veb,Nn and veb,En are the north and east components of the craft’s velocity, h is the craft’s altitude, and RE is earth’s radius.Furthermore,ϕi−1, θi−1, and ψi−1 are the previous epoch’s roll, pitch, and yaw sigma points and are the first three elements of each column of χ, and bp,q,r is the sigma points corresponding to the estimated IMU bias, which is predicted as random walk parameters.
(20)fbp,q,r:bp,ibq,ibr,i=bp,i−1bq,i−1br,i−1+wb,pwb,qwb,r.The measurement-prediction matrixΨ is populated by the predicted measurement vectors using each set of sigma-points in χ. Because measurements occur at different rates in this filter, it is necessary to have different measurement updates occur at different rates. For epochs coinciding with the sun sensor and GNSS attitude measurements, each column Ψi is as follows:
(21)Ψi=Bb,xBb,yBb,z∠X∠Yϕ′θ′ψ′,where Bb, ∠X, and ∠Y are predicted magnetometer and sun sensor measurements based on the ith sigma point.The observation equations useC^nb, the direction-cosine representation of the predicted attitude states ϕ^, θ^, and ψ^. Modeling of the magnetometer measurements uses the attitude state multiplied by the earth’s magnetic field vector:
(22)hBϕ,θ,ψ:B→b=C^nbB→n.Sun sensor measurements begin with the navigation-frame solar incidence vectorVsun,n, which are multiplied by the attitude state to find the body-frame incidence vector
(23)Vsun,b=C^nbVsun,n,and the solar incidence angles ∠X and ∠Y are then calculated as follows:
(24)h∠Xϕ,θ,ψ:∠X=π2+atan2Vsun,z,bVsun,x,b,h∠Yϕ,θ,ψ:∠Y=π2+atan2sunz,bsuny,b,where atan2 is the four-quadrant tangent inverse.The measurement update matrixz consists of the simulated sensor measurement at each filter epoch. These are similar in form to the columns of Ψ:
(25)z=Bb,xBb,yBb,z∠X∠YϕGNSSθGNSSψGNSS.For epochs without GNSS and sun sensor measurements,Ψ and z consist only of magnetometer predictions and measurements.


## 3. Data Simulation

### 3.1. Flight Profile

The simulated flight data used in this study is based upon the recorded flight data of ANITA III balloon experiment. That is, to simulate a balloon flight, the onboard position and attitude solutions were accepted as truth for simulation purposes, and sensor readings with realistic measurement noise were simulated.Figure3 shows the Euler angle time histories during a two-hour segment of the ANITA III flight. As indicated in Figure 3, the platform had a small (<1°) oscillation in the roll and pitch axes and a constant rotation about the yaw axis. A variable starting location was used to investigate the effect of the lower GDOP and ADOP at high latitudes.Figure 3
Attitude profile used in this work. This figure has been reproduced from the conference paper version of this work [16].


### 3.2. GNSS Observable Simulation

For each simulation run, four GNSS receivers were simulated with baseline separations of one-meter each, such that they are arranged in a square configuration. That is, the antennas were placed according to the following matrix:(26)Lb=x2,by2,bz2,bx3,by3,bz3,bx4,by4,bz4,b,where xi,b, yi,b, and zi,b are the body-centric coordinates of the ith antenna, i=1, denoting the master antenna. GNSS carrier-phase data was simulated for each flight profile at a rate of 10 Hz using the MATLAB SatNav toolbox [36], which was modified by Watson et al. [37] to include additional GNSS error sources.A number of deterministic and nondeterministic error sources are associated with GNSS measurements [21]. Fortunately, for attitude estimation applications, several of the primary GNSS error sources, including satellite and receiver clock biases and atmospheric delays, are canceled through the use of double-differenced GNSS observations [21]. However, two important error sources, namely, multipath reflections and carrier-phase breaks, or cycle slips, remain present. For example, when a metallic object reflects a GNSS signal onto the antenna, the multiple paths induce errors [21]. This could be a large problem on balloon-based scientific platforms, as the antennas are spaced closely and in close proximity to science payload. Thermal measurement noise in the receiver is another error source; it is actually amplified by double differencing GNSS data. As such, for this simulation study, multipath, carrier-phase breaks, and receiver thermal errors were assessed with respect to their effect on the attitude estimator’s performance using the distributions indicated in Table 2.Table 2
Sensor error source Monte Carlo simulation distribution parameters, same as in [16].
Error sources
Model parameters
Notes
Thermal noise
σ
ρ
=
0.32
m, σϕ=0.16λ
Linear scale factor randomly selected between 0 and 1
Multipath
1.0 intensity:σ=0.4m, τ=15sec
Linear scale factor randomly selected between 0 and 2
Tropospheric delay
Percent of error assumed handled by broadcast correction
Modified Hopfield with linear scale factor randomly selected between 0.95 and 1.05
Ionospheric delay
First-order ionospheric effects mitigated with dual frequency
Linear scale factor randomly selected between 0.7 and 1
Carrier phase break
Likelihood set to 1 phase break per 24 minutes to 1 phase break per 240 minutes
Gyroscope
In-run biasσ=9.6e−6rad/sec, ARW=0.2deg/hr
Scaled Honeywell HG1700AG72 SF = (1/50, 1/200, 1/400)
Sun sensor
Zenith measurement noiseσ=0.1deg
Scaled solar MEMS ISSDX-60 SF = (1, 2, 3, 4)
Magnetometer
Measurement noiseσ=2.67nTAsi terms scaled between 0.005 and 0.01B→hi terms scaled between 25 nT and 50 nT
ST LSM9DS0


### 3.3. Inertial Measurement Simulation

In addition to GNSS measurements, inertial measurement unit data was simulated for each flight profile and data at a sampling rate of 200 Hz. In particular, four grades of IMU triaxial rate gyroscope and accelerometers were assessed in the simulation. To simulate the effect of IMU noise, these ideal measurements were then polluted with both a time-varying biasbi with in-run stability σin‐run and measurement noise σARW:
(27)bi=bi−1+X1σin‐run,Δθi′=Δθi+bi+X2σARW,where X1 and X2 are normally distributed random numbers with a zero mean and unit standard deviation. The magnitude of these two noise terms was selected based on the grade of the inertial sensors assumed, which were varied as indicated in Table 2.


### 3.4. Magnetometer Simulation

Triaxial magnetometer data was also simulated for each flight based on the measurement models and uncertainties of available sensors [38]. In particular, the earth’s magnetic field along the flight profile was calculated and sensor measurements were simulated by polluting these true values with random noise based on the measurement noises quoted by the manufacturers’ specification sheets as indicated in Table 2 [38].The magnetometer data consists of magnetic field intensity measurements (Bb) in three orthogonal directions corresponding to the north, N, east, E, and down D axes in the body frame, b. This begins with BE, a vector constraining the simulated magnetic field intensities in the navigation frame, generated at each location along the flight path:
(28)B→E=Bb,NBb,EBb,D.The magnetic field was generated using the World Magnetic Model (WMM) [39] in an interface developed by Hardy et al. [40].Body-frame magnetic field measurements are simulated by multiplying truth attitude (represented by the direction-cosine matrixCnb) by the navigation-frame magnetic field:
(29)B→b=CnbB→E.With three contributing error sources added: hard and soft iron errors and measurement noise, in a simplified method as described by Gebre-Egziabher et al. [41]:
(30)B^=AsiB→b+B→hi,where Asi is a 3 × 3 matrix which describes the soft-iron error effect and B→hi is a 3 × 1 vector containing the hard-iron offset, a magnetic field generated by ferromagnetic material on the platform. For this study, nominal values for Asi and B→hi were used, based on the calibrations in the Gebre-Egziabher et al. paper [41]. Simulated measurement noise was then added to B^, corresponding to the precision level of the modeled magnetometer.


### 3.5. Sun Sensor Simulation

Simulated sun sensor data consists of solar incidence angles∠X and ∠Y relative to the two horizontal body-frame axes Xb and Yb. These were generated using the apparent solar azimuth θsun and elevation ϕsun calculated for each epoch of the flight duration. First, the solar azimuth and elevation values are transformed into a unit vector representing the sun’s position in the sky with respect to the navigation frame, n:
(31)Vsun,n=sunx,nsuny,nsunz,n.This unit-vector is then transformed using the nav-to-body direction cosine matrixCnb:
(32)Vsun,b=CnbVsun,n.The incidence angles∠X and ∠Y are calculated as follows:
(33)∠X=π2+atan2sunz,bsunx,b,∠Y=π2+atan2sunz,bsuny,b,as in the filter’s measurement prediction step.As with the magnetometer measurements, simulated measurement noise was added to the sun sensor measurements, corresponding to the performance of available sensors [42]. However, in the case of a sun sensor, as measurement noise increases at low solar elevations, the measurement noise was scaled according to solar elevation angle. Sun sensor measurements were simulated at 10 Hz intervals.


### 3.6. Error Sources and Monte Carlo

For this study, a total of 50 one-hour flight profiles were simulated in a Monte Carlo manner. In particular, the ECEF starting positions, magnitude of GNSS error sources, and quality of IMU, magnetometer, and sun sensor data were varied as indicated in Table2. Note that by randomly varying the starting location, the GNSS constellation satellite geometry was randomized as well.


## 3.1. Flight Profile

The simulated flight data used in this study is based upon the recorded flight data of ANITA III balloon experiment. That is, to simulate a balloon flight, the onboard position and attitude solutions were accepted as truth for simulation purposes, and sensor readings with realistic measurement noise were simulated.Figure3 shows the Euler angle time histories during a two-hour segment of the ANITA III flight. As indicated in Figure 3, the platform had a small (<1°) oscillation in the roll and pitch axes and a constant rotation about the yaw axis. A variable starting location was used to investigate the effect of the lower GDOP and ADOP at high latitudes.Figure 3
Attitude profile used in this work. This figure has been reproduced from the conference paper version of this work [16].


## 3.2. GNSS Observable Simulation

For each simulation run, four GNSS receivers were simulated with baseline separations of one-meter each, such that they are arranged in a square configuration. That is, the antennas were placed according to the following matrix:(26)Lb=x2,by2,bz2,bx3,by3,bz3,bx4,by4,bz4,b,where xi,b, yi,b, and zi,b are the body-centric coordinates of the ith antenna, i=1, denoting the master antenna. GNSS carrier-phase data was simulated for each flight profile at a rate of 10 Hz using the MATLAB SatNav toolbox [36], which was modified by Watson et al. [37] to include additional GNSS error sources.A number of deterministic and nondeterministic error sources are associated with GNSS measurements [21]. Fortunately, for attitude estimation applications, several of the primary GNSS error sources, including satellite and receiver clock biases and atmospheric delays, are canceled through the use of double-differenced GNSS observations [21]. However, two important error sources, namely, multipath reflections and carrier-phase breaks, or cycle slips, remain present. For example, when a metallic object reflects a GNSS signal onto the antenna, the multiple paths induce errors [21]. This could be a large problem on balloon-based scientific platforms, as the antennas are spaced closely and in close proximity to science payload. Thermal measurement noise in the receiver is another error source; it is actually amplified by double differencing GNSS data. As such, for this simulation study, multipath, carrier-phase breaks, and receiver thermal errors were assessed with respect to their effect on the attitude estimator’s performance using the distributions indicated in Table 2.Table 2
Sensor error source Monte Carlo simulation distribution parameters, same as in [16].
Error sources
Model parameters
Notes
Thermal noise
σ
ρ
=
0.32
m, σϕ=0.16λ
Linear scale factor randomly selected between 0 and 1
Multipath
1.0 intensity:σ=0.4m, τ=15sec
Linear scale factor randomly selected between 0 and 2
Tropospheric delay
Percent of error assumed handled by broadcast correction
Modified Hopfield with linear scale factor randomly selected between 0.95 and 1.05
Ionospheric delay
First-order ionospheric effects mitigated with dual frequency
Linear scale factor randomly selected between 0.7 and 1
Carrier phase break
Likelihood set to 1 phase break per 24 minutes to 1 phase break per 240 minutes
Gyroscope
In-run biasσ=9.6e−6rad/sec, ARW=0.2deg/hr
Scaled Honeywell HG1700AG72 SF = (1/50, 1/200, 1/400)
Sun sensor
Zenith measurement noiseσ=0.1deg
Scaled solar MEMS ISSDX-60 SF = (1, 2, 3, 4)
Magnetometer
Measurement noiseσ=2.67nTAsi terms scaled between 0.005 and 0.01B→hi terms scaled between 25 nT and 50 nT
ST LSM9DS0


## 3.3. Inertial Measurement Simulation

In addition to GNSS measurements, inertial measurement unit data was simulated for each flight profile and data at a sampling rate of 200 Hz. In particular, four grades of IMU triaxial rate gyroscope and accelerometers were assessed in the simulation. To simulate the effect of IMU noise, these ideal measurements were then polluted with both a time-varying biasbi with in-run stability σin‐run and measurement noise σARW:
(27)bi=bi−1+X1σin‐run,Δθi′=Δθi+bi+X2σARW,where X1 and X2 are normally distributed random numbers with a zero mean and unit standard deviation. The magnitude of these two noise terms was selected based on the grade of the inertial sensors assumed, which were varied as indicated in Table 2.


## 3.4. Magnetometer Simulation

Triaxial magnetometer data was also simulated for each flight based on the measurement models and uncertainties of available sensors [38]. In particular, the earth’s magnetic field along the flight profile was calculated and sensor measurements were simulated by polluting these true values with random noise based on the measurement noises quoted by the manufacturers’ specification sheets as indicated in Table 2 [38].The magnetometer data consists of magnetic field intensity measurements (Bb) in three orthogonal directions corresponding to the north, N, east, E, and down D axes in the body frame, b. This begins with BE, a vector constraining the simulated magnetic field intensities in the navigation frame, generated at each location along the flight path:
(28)B→E=Bb,NBb,EBb,D.The magnetic field was generated using the World Magnetic Model (WMM) [39] in an interface developed by Hardy et al. [40].Body-frame magnetic field measurements are simulated by multiplying truth attitude (represented by the direction-cosine matrixCnb) by the navigation-frame magnetic field:
(29)B→b=CnbB→E.With three contributing error sources added: hard and soft iron errors and measurement noise, in a simplified method as described by Gebre-Egziabher et al. [41]:
(30)B^=AsiB→b+B→hi,where Asi is a 3 × 3 matrix which describes the soft-iron error effect and B→hi is a 3 × 1 vector containing the hard-iron offset, a magnetic field generated by ferromagnetic material on the platform. For this study, nominal values for Asi and B→hi were used, based on the calibrations in the Gebre-Egziabher et al. paper [41]. Simulated measurement noise was then added to B^, corresponding to the precision level of the modeled magnetometer.


## 3.5. Sun Sensor Simulation

Simulated sun sensor data consists of solar incidence angles∠X and ∠Y relative to the two horizontal body-frame axes Xb and Yb. These were generated using the apparent solar azimuth θsun and elevation ϕsun calculated for each epoch of the flight duration. First, the solar azimuth and elevation values are transformed into a unit vector representing the sun’s position in the sky with respect to the navigation frame, n:
(31)Vsun,n=sunx,nsuny,nsunz,n.This unit-vector is then transformed using the nav-to-body direction cosine matrixCnb:
(32)Vsun,b=CnbVsun,n.The incidence angles∠X and ∠Y are calculated as follows:
(33)∠X=π2+atan2sunz,bsunx,b,∠Y=π2+atan2sunz,bsuny,b,as in the filter’s measurement prediction step.As with the magnetometer measurements, simulated measurement noise was added to the sun sensor measurements, corresponding to the performance of available sensors [42]. However, in the case of a sun sensor, as measurement noise increases at low solar elevations, the measurement noise was scaled according to solar elevation angle. Sun sensor measurements were simulated at 10 Hz intervals.


## 3.6. Error Sources and Monte Carlo

For this study, a total of 50 one-hour flight profiles were simulated in a Monte Carlo manner. In particular, the ECEF starting positions, magnitude of GNSS error sources, and quality of IMU, magnetometer, and sun sensor data were varied as indicated in Table2. Note that by randomly varying the starting location, the GNSS constellation satellite geometry was randomized as well.


## 4. Results

### 4.1. GNSS Only

The GNSS-only attitude determination script was run in two modes, the first using GPS data only and the second adding GLONASS observables. The pitch, roll, and heading error statistics for both filter modes are presented in Table3. These results include two simulations for which the baseline filter solution failed to converge, presumably due to carrier-phase break.Table 3
GNSS attitude performance—median RMS error.
Median errors
Roll (deg.)
Pitch (deg.)
Heading (deg.)
GPS only
0.44
0.42
0.16
GPS + GLONASS
0.28
0.26
0.12Using GLONASS as well as GPS satellites yielded a median performance improvement of 40 percent lower attitude error. In an Antarctic flight regime, fewer GNSS satellites are observable, and these are seen at lower elevations [43]. This can negatively impact the geometric dilution of precision (GDOP), a metric that describes the geometric diversity of satellite-receiver vectors [21] and also the attitude dilution of precision (ADOP), as defined in (14). Figure 4 shows the ADOP calculated by latitude for a single meridian, averaged over a full day. At high latitudes, the GLONASS constellation has a lower (and therefore more desirable) ADOP figure, and using both constellations yields a lower ADOP at all latitudes. Figure 5 shows the overall comparison between the GPS-only and GPS and GLONASS attitude solutions. For most profiles across the latitude range, the GPS + GLONASS solution outperformed the GPS-only solution.Figure 4
ADOP calculated along one meridian, averaged over a 24-hour period for the GPS constellation, for the GLONASS constellation, and for both constellations.Figure 5
Comparison between GPS-only mode and GPS + GLONASS mode for all profiles, latitude shown.


### 4.2. GNSS and Multisensor Attitude Filter

Table4 presents overall error statistics for the 50 trials for the GNSS + INS, GNSS + all sensors, and all sensors without GNSS, respectively.Table 4
Unscented Kalman filter error statistics: median attitude error.
Median errors
Roll (deg.)
Pitch (deg.)
Heading (deg.)
GNSS + INS
0.59
0.38
0.64
INS + Mag + SS
0.0041
0.047
0.045
GNSS + INS + Mag + SS
0.048
0.053
0.047Figure6 shows the cumulative distribution of the 3D attitude error = ϕ2+θ2+ψ2 for the various filter configurations over the 50 simulated flights, and Figure 7 shows the corresponding roll, pitch, and yaw errors for the simulated flights.Figure 6
Comparison between GNSS-SVD solution and multisensor attitude filter in different modes—3-axis attitude error.Figure 7
Comparison between GNSS-SVD solution and multisensor attitude filter in different modes—roll, pitch, and yaw error.In these tables, it is clear that using additional sensors in addition to GNSS can markedly improve performance. For example, Figure8 shows the attitude estimation error for one example trial, in which the GNSS-only attitude is shown alongside the multisensor filters for comparison. When accurate pitch and roll estimates are critical, the inertial sensors offer more trials with significant accuracy increases, as evident when looking at the GNSS-INS-only solution. When the full suite of sensors is combined, 10% of the sun/INS/Mag solutions that drift considerably are improved. As expected, the yaw estimate is greatly improved by the sun and magnetic measurements.Figure 8
Roll, pitch, and heading errors for multisensor filter in GNSS + INS mode, GNSS + INS + Mag + SS mode, with GNSS-only result for comparison. This figure has been reproduced from the conference paper version of this work [16].Of great interest is the algorithm’s ability to handle carrier-phase breaks. For example, phase breaks could occur due to radiofrequency interference, such as during a data uplink/downlink transmission over the Iridium satellite constellation which operates very close to the GPS L1 frequency [44]. When a carrier-phase break occurs, it can fortunately be detected easily by a data editor [45]. As such, whenever this occurs, the baseline estimation filter resets the error covariance for the impacted carrier-phase ambiguities to a large value. The multisensor filter attitude determination performance was generally lower across the range of phase break likelihoods as shown in Figure 9. Notably, the multisensor UKF yielded a low error-level attitude solution for the few trials with considerably high GNSS-attitude errors.Figure 9
RMS attitude versus phase break likelihood for each trial.Also of interest is the estimator’s performance with high receiver measurement thermal noise and GNSS multipath reflection errors. These could arise on a scientific platform depending on the GNSS antenna configuration relative to the other scientific instruments and antennas/Figures10 and 11 show that the multisensor filter yields lower magnitude errors than the GNSS-only filter across both error scale ranges. Although an increasing level of multipath error did not noticeably affect the result of the GNSS-only filter performance, the multisensor filter performed better in nearly all trials.Figure 10
RMS attitude versus thermal error scale factor for each trial.Figure 11
RMS attitude versus multipath error scale factor for each trial.Sensitivity to the ionospheric and tropospheric error contribution to the GNSS errors was not considered, as the short baseline between antennas led to cancellation of those error sources.


## 4.1. GNSS Only

The GNSS-only attitude determination script was run in two modes, the first using GPS data only and the second adding GLONASS observables. The pitch, roll, and heading error statistics for both filter modes are presented in Table3. These results include two simulations for which the baseline filter solution failed to converge, presumably due to carrier-phase break.Table 3
GNSS attitude performance—median RMS error.
Median errors
Roll (deg.)
Pitch (deg.)
Heading (deg.)
GPS only
0.44
0.42
0.16
GPS + GLONASS
0.28
0.26
0.12Using GLONASS as well as GPS satellites yielded a median performance improvement of 40 percent lower attitude error. In an Antarctic flight regime, fewer GNSS satellites are observable, and these are seen at lower elevations [43]. This can negatively impact the geometric dilution of precision (GDOP), a metric that describes the geometric diversity of satellite-receiver vectors [21] and also the attitude dilution of precision (ADOP), as defined in (14). Figure 4 shows the ADOP calculated by latitude for a single meridian, averaged over a full day. At high latitudes, the GLONASS constellation has a lower (and therefore more desirable) ADOP figure, and using both constellations yields a lower ADOP at all latitudes. Figure 5 shows the overall comparison between the GPS-only and GPS and GLONASS attitude solutions. For most profiles across the latitude range, the GPS + GLONASS solution outperformed the GPS-only solution.Figure 4
ADOP calculated along one meridian, averaged over a 24-hour period for the GPS constellation, for the GLONASS constellation, and for both constellations.Figure 5
Comparison between GPS-only mode and GPS + GLONASS mode for all profiles, latitude shown.


## 4.2. GNSS and Multisensor Attitude Filter

Table4 presents overall error statistics for the 50 trials for the GNSS + INS, GNSS + all sensors, and all sensors without GNSS, respectively.Table 4
Unscented Kalman filter error statistics: median attitude error.
Median errors
Roll (deg.)
Pitch (deg.)
Heading (deg.)
GNSS + INS
0.59
0.38
0.64
INS + Mag + SS
0.0041
0.047
0.045
GNSS + INS + Mag + SS
0.048
0.053
0.047Figure6 shows the cumulative distribution of the 3D attitude error = ϕ2+θ2+ψ2 for the various filter configurations over the 50 simulated flights, and Figure 7 shows the corresponding roll, pitch, and yaw errors for the simulated flights.Figure 6
Comparison between GNSS-SVD solution and multisensor attitude filter in different modes—3-axis attitude error.Figure 7
Comparison between GNSS-SVD solution and multisensor attitude filter in different modes—roll, pitch, and yaw error.In these tables, it is clear that using additional sensors in addition to GNSS can markedly improve performance. For example, Figure8 shows the attitude estimation error for one example trial, in which the GNSS-only attitude is shown alongside the multisensor filters for comparison. When accurate pitch and roll estimates are critical, the inertial sensors offer more trials with significant accuracy increases, as evident when looking at the GNSS-INS-only solution. When the full suite of sensors is combined, 10% of the sun/INS/Mag solutions that drift considerably are improved. As expected, the yaw estimate is greatly improved by the sun and magnetic measurements.Figure 8
Roll, pitch, and heading errors for multisensor filter in GNSS + INS mode, GNSS + INS + Mag + SS mode, with GNSS-only result for comparison. This figure has been reproduced from the conference paper version of this work [16].Of great interest is the algorithm’s ability to handle carrier-phase breaks. For example, phase breaks could occur due to radiofrequency interference, such as during a data uplink/downlink transmission over the Iridium satellite constellation which operates very close to the GPS L1 frequency [44]. When a carrier-phase break occurs, it can fortunately be detected easily by a data editor [45]. As such, whenever this occurs, the baseline estimation filter resets the error covariance for the impacted carrier-phase ambiguities to a large value. The multisensor filter attitude determination performance was generally lower across the range of phase break likelihoods as shown in Figure 9. Notably, the multisensor UKF yielded a low error-level attitude solution for the few trials with considerably high GNSS-attitude errors.Figure 9
RMS attitude versus phase break likelihood for each trial.Also of interest is the estimator’s performance with high receiver measurement thermal noise and GNSS multipath reflection errors. These could arise on a scientific platform depending on the GNSS antenna configuration relative to the other scientific instruments and antennas/Figures10 and 11 show that the multisensor filter yields lower magnitude errors than the GNSS-only filter across both error scale ranges. Although an increasing level of multipath error did not noticeably affect the result of the GNSS-only filter performance, the multisensor filter performed better in nearly all trials.Figure 10
RMS attitude versus thermal error scale factor for each trial.Figure 11
RMS attitude versus multipath error scale factor for each trial.Sensitivity to the ionospheric and tropospheric error contribution to the GNSS errors was not considered, as the short baseline between antennas led to cancellation of those error sources.


## 5. Conclusions

This study outlined the design of a GNSS-based attitude determination algorithm and its incorporation in a nonlinear attitude-determination filter which also utilized inertial measurements. Additional measurements from sun sensors and magnetometers were also proposed and added to the algorithm architecture. The algorithm was run in multiple modes with varying levels of measurement errors to assess how each measurement type contributed to the attitude determination performance. First, it was shown that the GNSS-only attitude solutions are consistently improved when GLONASS satellites are included in addition to GPS, owing to more observables and lower dilution of precision (especially in polar regions). Across all flight profiles, GLONASS observables yielded a median of 30 percent lower error in the pitch and roll axes, with closer performance in the yaw axis. The performance was somewhat mitigated when GNSS Euler angle measurements were combined with inertial measurements only. This was most evident in the yaw axis, further indicating tuning as a potential culprit. Incorporating sun sensor and magnetometer measurements yielded the best improvement to the nonlinear filter’s performance, with median performance better than tenth-degree accuracy for a majority of the trials. The combined estimator showed nearly a uniform distribution and nearly consistent improvement when assessing sensitivity against two important GNSS error sources, phase breaks and multipath.



---

*Source: 4871239-2018-01-16.xml*