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.
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 Section 2, 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
Figure 1 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.
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, , for this filter consists of the relative position vector components between antennas and , , , and , and a set of double-differenced pseudoranges and carrier-phase biases, .
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] where is the wavelength corresponding to the frequencies and and expressed in meters. The geometric range between the receiver and GNSS satellite is also expressed in meters, as are the ionospheric and tropospheric delays and . The speed of light is expressed in meters per second. The clock biases of the receiver and satellite, and , 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 antenna (antenna 1) and (antennas 2, 3, or 4) are differenced to form single-differenced phase measurements:
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 and reference satellite , where the remaining receiver clock bias errors are eliminated, leaving only the unknown phase biases , 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, , transforms the state 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, .
The measurement vector, , consists of double-differenced pseudorange and phase measurements for each satellite relative to the reference satellite, including measurements for each and frequencies:
Operating in parallel with the baseline estimation Kalman filter, the floating point estimated phase biases (for GPS satellites only); 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, , and measurement noise, , and initial error-covariance, , assumptions selected for the differential GNSS baseline estimator are outlined in Table 1:
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, , is constructed at each epoch by vertically concatenating the estimate relative vectors of each of nonmaster antenna, as adopted from Cohen and Parkinson [3]:
This matrix is used to estimate the platform attitude given by the antenna baseline vectors. The known bodycentric antenna coordinate matrix, , 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 between the ECEF, , and body, , frames [26]. This is then transformed into the local navigation frame using, , the earth-to-nav rotation, which is dependent on the master antennas’ latitude and longitude and is defined in many texts [20]:
This solution requires the construction of a matrix using the measured vectors and reference vectors :
A SVD is performed on , resulting in unitary matrices and . A diagonal, 3 × 3 matrix is constructed using the determinants of and :
By multiplying , , and , one can find the rotation matrix :
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 matrix is multiplied by the transpose of the nontransformed rotation matrix .
can then be used to find the inverse of the error covariance matrix:
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] where is the number of satellites in view, is the 3 × 3 identity matrix, and 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 , and observation functions 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, estimated in the multisensor filter is given as where , , and are the platform’s roll, pitch, and yaw and is the time-varying biases of the IMU’s roll rate, , pitch rate, , and yaw rate, , gyroscopes.
Within the UKF framework, at each epoch, the state vector is expanded into a group of sigma points, , where is the length of the estimated state vector. For each group of sigma points , the attitude states are predicted by integrating the IMU gyro data through the attitude kinematic equations [35]: where represents sine, represents cosine, and represents tangent.
The delta angles , , and , 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.
The earth-rate rotations , , and can be found using the earth’s rotation transformed to the body frame, then converting to Euler angles [20] where is earth’s rotation rate and is the craft’s latitude. The craft-rate rotation component can be found in a similar way [20] where and are the north and east components of the craft’s velocity, is the craft’s altitude, and is earth’s radius.
Furthermore, , , and are the previous epoch’s roll, pitch, and yaw sigma points and are the first three elements of each column of , and is the sigma points corresponding to the estimated IMU bias, which is predicted as random walk parameters.
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 is as follows: where , , and are predicted magnetometer and sun sensor measurements based on the ith sigma point.
The observation equations use , 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:
Sun sensor measurements begin with the navigation-frame solar incidence vector , which are multiplied by the attitude state to find the body-frame incidence vector and the solar incidence angles and are then calculated as follows: where is the four-quadrant tangent inverse.
The measurement update matrix consists of the simulated sensor measurement at each filter epoch. These are similar in form to the columns of :
For epochs without GNSS and sun sensor measurements, and 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.
Figure 3 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.
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: where , , and are the body-centric coordinates of the ith antenna, , 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.
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 bias with in-run stability and measurement noise : where and 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 () in three orthogonal directions corresponding to the north, , east, , and down axes in the body frame, . This begins with , a vector constraining the simulated magnetic field intensities in the navigation frame, generated at each location along the flight path:
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 matrix ) by the navigation-frame magnetic field:
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]: where is a 3 × 3 matrix which describes the soft-iron error effect and 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 and were used, based on the calibrations in the Gebre-Egziabher et al. paper [41]. Simulated measurement noise was then added to , corresponding to the precision level of the modeled magnetometer.
3.5. Sun Sensor Simulation
Simulated sun sensor data consists of solar incidence angles and relative to the two horizontal body-frame axes and . These were generated using the apparent solar azimuth and elevation 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, :
This unit-vector is then transformed using the nav-to-body direction cosine matrix :
The incidence angles and are calculated as follows: 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 Table 2. 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 Table 3. These results include two simulations for which the baseline filter solution failed to converge, presumably due to carrier-phase break.
Using 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.
4.2. GNSS and Multisensor Attitude Filter
Table 4 presents overall error statistics for the 50 trials for the GNSS + INS, GNSS + all sensors, and all sensors without GNSS, respectively.
Figure 6 shows the cumulative distribution of the 3D attitude error = 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.
In these tables, it is clear that using additional sensors in addition to GNSS can markedly improve performance. For example, Figure 8 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.
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.
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/Figures 10 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.
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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by NASA West Virginia Space Grant Consortium Graduate Student Fellowship Program and the National Geospatial-Intelligence Agency Academic Research Program Grant no. HM0476-15-1-0004. This study is approved for NGA public release under case no. 17-839.