Calibration of Magnetometers with GNSS Receivers and Magnetometer-Aided GNSS Ambiguity Fixing

Magnetometers provide compass information, and are widely used for navigation, orientation and alignment of objects. As magnetometers are affected by sensor biases and eventually by systematic distortions of the Earth magnetic field, a calibration is needed. In this paper, a method for calibration of magnetometers with three Global Navigation Satellite System (GNSS) receivers is presented. We perform a least-squares estimation of the magnetic flux and sensor biases using GNSS-based attitude information. The attitude is obtained from the relative positions between the GNSS receivers in the North-East-Down coordinate frame and prior knowledge of these relative positions in the platform’s coordinate frame. The relative positions and integer ambiguities of the periodic carrier phase measurements are determined with an integer least-squares estimation using an integer decorrelation and sequential tree search. Prior knowledge on the relative positions is used to increase the success rate of ambiguity fixing. We have validated the proposed method with low-cost magnetometers and GNSS receivers on a vehicle in a test drive. The calibration enabled a consistent heading determination with an accuracy of five degrees. This precise magnetometer-based attitude information allows an instantaneous GNSS integer ambiguity fixing.


Introduction
A precise attitude information is essential for numerous applications, e.g., for navigation, for flight control and for constructions. Magnetometers provide attitude information globally but are affected by sensor biases and measurement errors. The latter ones can be systematic due to metallic environment and/or stochastic due to measurement noise. Sensor biases and systematic errors can be removed by calibration.
The calibration can be performed with three Global Navigation Satellite System (GNSS) receivers (e.g., GPS, Galileo, GLONASS, Beidou) being placed on the same platform as the magnetometer. Satellite navigation enables the determination of the relative position between both GNSS receivers in a local coordinate frame, e.g., the North-East-Down (NED) frame. The obtained relative position is related to the prior knowledge of the relative position in the body-fixed frame to obtain the platform's attitude.
The use of GNSS carrier phase measurements is required to obtain a precise attitude information. The carrier phases are periodic, which requires an integer ambiguity resolution for each differential measurement. An efficient integer least-squares estimation of these integer ambiguities has been developed by Teunissen in [1]. His "Least Squares Ambiguity Decorrelation Adjustment" (LAMBDA) method includes an integer decorrelation and sequential tree search. The success rate of integer ambiguity fixing can be substantially increased by including prior information on the baseline. Teunissen has included prior information on the baseline length in [2,3] and Henkel and Günther have additionally included prior information on the attitude in [4]. The prior information on the length and attitude angles can be either a "hard" or "soft" information. Henkel et al. have developed a method for determining magnetometer biases using two GNSS receivers in [5]. Crassidis et al. [6] and Psiaki et al. [7] used extended and unscented Kalman filters to estimate the sensor biases and scaling factors.
The purpose of this paper is threefold: first, we describe a method for calibration of magnetometer biases using attitude information from three GNSS receivers. Second, we describe a method for fast and reliable integer ambiguity fixing using magnetometer-based attitude information. Finally, the proposed methods are verified with measurements of low-cost magnetometers and GNSS receivers in a test drive with a vehicle.

Measurement Models
In this section, we introduce a precise model for the magnetometer measurements.

Measurement Model for Magnetic Field Sensor
We model the magnetic flux measurement in the sensor-fixed s-frame at epoch t n as a sum of true magnetic flux, measurement bias b m s and measurement noise ε m s : where the magnetic flux is expressed in the local geomagnetic (g-) frame (right-hand coordinate system with the first axis pointing towards geomagnetic north pole, the second axis parallel to plane of magnetic flux, and the third axis perpendicular to magnetic flux and pointing downwards) and transformed into the sensor-fixed frame. The transformation R s g is performed in three steps: • transformation from local geomagnetic frame into local geographic navigation-frame, depending on the magnetic declination δν, • transformation from local geographic navigation-frame into body-fixed frame, depending on the attitude of sensor's platform (roll ϕ, pitch θ and heading ψ), • transformation from body-fixed frame into sensor-fixed frame, depending on misalignment errors of sensors (roll offset ∆ϕ, pitch offset ∆θ and heading offset ∆ψ).

Calibration of Magnetic Flux Sensors
The calibration of the magnetic flux sensor requires the estimation of the magnitude m g of the magnetic flux and of the sensor bias b m s . We consider the rotation matrix R s g as known. It is derived from GNSS in this section.

Estimation of Magnetometer Biases
The calibration of the magnetic flux requires measurements with rotational dynamics of multiple epochs to enable a separation of m g and b m s . The magnetic flux m g and bias b m s can be assumed constant over a few epochs. Therefore, we stack the magnetic flux measurement of Equation (1) of n epochs in a column vector and express them in terms of the unknowns: The magnetic flux and bias are determined such that the sum of squared measurement residuals is minimized, i.e.,

Attitude Determination with Three GNSS Receivers
In this section, we derive the attitude information needed for calibration of the magnetometer from three GNSS receivers.
It is assumed that all of the three GNSS receivers and the magnetometer are mounted on the same platform. The first GNSS receiver serves as reference receiver for attitude determination. The reference receiver and the other two receivers span two independent baselines pointing from the reference receiver to the other receivers. GNSS enables the estimation of these baselines in a local North-East-Down coordinate frame. The relation of the baseline in the North-East-Down and local body-fixed coordinate frames provides the attitude information.
GNSS receivers provide two types of range measurements with the following characteristics: • Carrier phase measurement λϕ k r , ⊕ carrier phase can be tracked with millimeter accuracy, carrier phase is period with λ = 19 cm and requires ambiguity resolution, • Pseudorange measurement ρ k r , ⊕ pseudorange is an unambiguous range measurement, pseudorange measurement is more sensitive to multipath, pseudorange measurement can only be tracked with meter-level accuracy.

Modeling of Differential GNSS Measurements
We use both types of measurements and perform differential measurements between the reference receiver (being indexed by 1) and any other receiver (being indexed by r) to eliminate atmospheric delays, orbital errors, satellite clock offsets and biases. Additionally, a reference satellite l is selected and the differential measurements of the reference satellite are subtracted from the differential measurements of any other satellite k ∈ {1, . . . , K} to eliminate receiver clock errors and biases. The obtained double difference carrier phase measurement is modeled as with the wavelength λ, the carrier phase measurement ϕ k r in unit of cycles, the normalized line of sight vector e k from satellite k to the receiver platform, the rotation matrix R e n from the navigation frame into the Earth Centered Earth Fixed (ECEF) (e-) coordinate frame depending on the absolute position of the platform (given by longitude λ 1 and latitude ϕ 1 ), the integer ambiguity N k r , the phase noise ε k r and the double difference operator (·) kl 1r = (·) k 1 − (·) k r − (·) l 1 − (· r ) l . Similarly, the double difference pseudorange measurement is modeled as with the pseudorange multipath error ∆ρ k MP,r and pseudorange noise η k r .

Joint Estimation of Baselines, Pseudorange Multipaths and Ambiguities
The double difference carrier phase and pseudorange measurements from both baselines are stacked in a common measurement vector: with The number of unknowns of Equation (10) exceeds the number of measurements of one epoch, but the ambiguities are constant over time. Thus, measurements from multiple epochs are required to estimate the baselines, ambiguities and pseudorange multipaths. We introduce a state space/movement model to describe the temporal behavior of the state parameters, i.e., with state transition matrix Φ n and process noise η x n . A Kalman filter (see Brown and Hwang [8]) is used to estimate the state vector. We first disregard the integer property of ambiguities to obtain a float solution. The Kalman filter includes the state prediction Sensors 2017, 17, 1324 5 of 13 and the subsequent state updatê

Integer Ambiguity Fixing Using Prior Information on Baseline Coordinates
In this section, we exploit the integer property of ambiguities and describe the fixing of float ambiguities to integer ones. A joint fixing of ambiguities from both baselines is performed to exploit the correlation introduced by the use of a common receiver in both baselines. The integer ambiguities, baselines and pseudorange multipaths are determined such that the sum of squared residuals is minimized, i.e., An analytical optimization is not feasible due to the multi-dimensional integer parameters. Therefore, a sequential tree search is performed as suggested by Teunissen in [1] and Jonge and de Tiberius [9]. The ambiguity residuals r N kl 1r :=N kl 1r − N kl 1r are expressed in terms of conditional ambiguitiesN kl|1l,...,(k−1)l 1r using the triangular decomposition of the float ambiguity covariance matrix (see Blewitt [10]). The lower and upper bound of the search interval of ambiguity N kl 1r are derived by Jonge and de Tiberius in [9] as: with χ 2 being the search space volume. A prior knowledge on the baselines is typically available in the body-fixed (b-) frame. The local (East-North-Up) navigation (n-) frame is related to the body-fixed frame via the attitude angles (roll ϕ, pitch θ and heading ψ), i.e., the baseline can be re-parameterized as The estimation of all three attitude angles requires at least three receivers spanning two baselines b 12 and b 13 . Once baseline estimates are available in both frames, the rotation matrix R n b (ϕ, θ, ψ) can be determined using Wahba's solution [11]. It minimizes the cost function Once the rotation matrix is found, the attitude angles can obtained from Equation (4): We constrain the sequential tree search of inequalitŷ by integrating prior information on the baseline lengths and coordinates (provided in body-fixed frame), i.e., we consider only these integer candidates N kl 1r , where the respective baseline estimate shows consistent baseline length and attitude information. This leads to the following additional requirements: with the following notations: We derive the baseline estimate for partially fixed ambiguities: the double difference ambiguities are split into float ambiguities N 1 1r and fixed ambiguities N 2 1r . As float and fixed double difference ambiguities in N 1r are not sorted, the mapping of N 1 1r and N 2 1r to N 1r is given by where A 1 1r and A 2 1r describe the mapping of float/integer ambiguity vectors N 1 1r and N 2 1r to N 1r . The ambiguity decomposition is included in the measurement model of Equation (10) being rewritten as The fixed ambiguitiesŇ 2 1r , r ∈ {2, 3} are subtracted from the measurements. Subsequently, a projection on the space orthogonal toĀ 1 is applied to eliminate all non-baseline terms, i.e., Sensors 2017, 17, 1324 7 of 13 The notation is further simplified by definingĀ := P ⊥ A 1H . The least-squares estimate of both partially fixed baselines follows as The accuracy of the baseline estimates depends on the number of fixed ambiguities: if less than three ambiguities are fixed per baseline, the pseudoranges are needed and the accuracy is limited by the pseudorange multipath. If at least three ambiguities are fixed per baseline, the baselines can be estimated solely with carrier phase measurements resulting in a much higher accuracy.
The efficiency of the ambiguity fixing is substantially improved by the integer decorrelation of Teunissen [1]. The integer decorrelation transformation Z is applied to the double difference ambiguities, i.e., The integer decorrelation transformation is invertible and preserves the integer property, i.e., the uncorrelated double difference ambiguities can be obtained from the decorrelated ones by The decorrelated fixed ambiguitiesÑ 1 1r and float ambiguitiesÑ 2 1r are introduced in the measurement model of Equation (25), which is rewritten as The decorrelated float ambiguities are eliminated by a projection on the space orthogonal ofÃ. The least-squares baseline estimates for partially fixed decorrelated ambiguities follows as withĀ = P ⊥ AH . The baseline estimates are used to constrain the search interval in Equations (22) and (23). The reduction of search space enables a much more reliable ambiguity fixing.

Fast Initialization of GNSS Attitude Ambiguity Fixing with Calibrated Magnetometers
In this section, we describe the benefit of magnetometer-based attitude information for kinematic GNSS-based attitude determination.
In a first step, we derive the attitude from calibrated magnetic flux measurements. It is assumed that precise estimates of the sensor biasb m s and of the magnetic flux strengthm g are available, e.g., using a previous calibration. In this case, the attitude information is obtained from the magnetic flux measurements by minimizing the sum of squared residuals, i.e., This minimization refers to the well-known Wahba's problem. Its solution is provided in [11]. In a second step, we use this attitude information as prior knowledge in the GNSS-based attitude determination. The GNSS carrier phase ambiguities are again determined with a sequential tree search as described in the previous section.
However, the second constraint (23) on the attitude substantially simplifies as the attitude angles are known from the calibrated magnetometer. Thus, the minimization and Wahba's solution [11] are no longer needed and Equation (23) simplifies to Figure 1 visualizes the constrained tree search. The integer candidates are determined sequentially for the satellites, whereas the integer candidates of each satellite are conditioned on the integer ambiguity candidates of the previous satellites. There are typically multiple integer candidates for each satellite. However, several integer candidates can be dropped due to inconsistent baseline length and attitude information, which reduces the number of branches to be searched.

Analysis of Benefit of Magnetometer-Based Attitude Information for GNSS Integer Ambiguity Fixing
In this section, we analyze the benefit of magnetometer-based attitude information for GNSS integer ambiguity fixing. Galileo double difference carrier phase and pseudorange measurements are simulated according to Equations (8) and (9) with the settings of Table 1.  Integer ambiguity fixing with sequential tree search and integer decorrelation using magnetometer-based attitude information and baseline length prior information. Figure 2 shows the probability of incorrect integer ambiguity fixing as a function of the accuracy of the magnetometer-based heading and pitch angle accuracies. We can observe that the probability of incorrect integer ambiguity fixing is reduced by several orders of magnitude by the magnetometer based heading/pitch angle with an accuracy of 10 • .

Measurement Results
This section describes the verification of the proposed magnetometer calibration and the achievable heading accuracy in a test drive using a tightly coupled GNSS/INS reference system. The method and results differ in the following aspects from our earlier paper [5]: The MPU 9250 [13] includes a three-axis silicon monolithic Hall-effect magnetic sensor with magnetic concentrator. The magnetic flux measurements are provided with a resolution of 0.6 µT. Three multi-sensor modules were mounted on the roof of a vehicle, whereas each multi-sensor module included the above sensors. The inertial sensor provides 3D acceleration, 3D angular rate and 3D magnetometer measurements in the sensor-fixed frame with an update rate of 100 Hz. The distances between the GNSS antennas were around 1 m. The Multi-GNSS (GPS + GLONASS)/INS tightly coupled attitude determination included some pre-processing (synchronization, cycle slip correction [14,15] calibration of inertial sensor biases) and carrier phase integer ambiguity fixing. Figure 4 shows the track of the test drive. It was split into two phases, i.e., a first phase for calibrating the magnetometer, and a second phase for testing the calibration of the magnetometer. In both phases, the track includes rotational dynamics. Figure 5 shows the accuracy of the magnetic flux measurements and of our measurement model, i.e., the raw measurements are compared to estimated measurements derived from the least-squares estimate of the magnetic flux in the North-East-Down frame and the sensor biases. Obviously, both tracks fit quite well. The systematic variations over time in both subfigures are caused by changes in the attitude (especially heading). The estimated measurements are less noisy as the precise Multi-GNSS/INS tightly coupled attitude is used and as constant magnitudes of the magnetic field in North-East-Down directions are assumed over the considered period of 40 s.          Figure 7. Detailed analysis of heading performance for two sections with moderate to high rotational dynamics: the filtered magnetometer-based heading deviates by less than 10 degrees from the GNSS/INS tightly coupled heading.

Conclusions
In this paper, a method for calibration of magnetic field sensors with three Global Navigation Satellite System (GNSS) receivers was provided. A GNSS-based attitude information was obtained from the relative positions between the GNSS receivers. The relative positions were jointly estimated with the carrier phase integer ambiguities and pseudorange multipath errors. The latter ones were considered as state parameters to exploit the spatial and temporal correlations of multipath. The reliability of carrier phase integer ambiguity fixing was improved by integration of prior information on the baseline coordinates provided in a local body-frame.
The attitude information of the calibrated magnetometer was also integrated into the sequential tree search of ambiguity fixing, i.e., integer candidates with inconsistent attitude information were Figure 7. Detailed analysis of heading performance for two sections with moderate to high rotational dynamics: the filtered magnetometer-based heading deviates by less than 10 degrees from the GNSS/INS tightly coupled heading.

Conclusions
In this paper, a method for calibration of magnetic field sensors with three Global Navigation Satellite System (GNSS) receivers was provided. A GNSS-based attitude information was obtained from the relative positions between the GNSS receivers. The relative positions were jointly estimated with the carrier phase integer ambiguities and pseudorange multipath errors. The latter ones were considered as state parameters to exploit the spatial and temporal correlations of multipath. The reliability of carrier phase integer ambiguity fixing was improved by integration of prior information on the baseline coordinates provided in a local body-frame.
The attitude information of the calibrated magnetometer was also integrated into the sequential tree search of ambiguity fixing, i.e., integer candidates with inconsistent attitude information were disregarded. It was shown that a heading and pitch angle accuracy of 10 degrees are sufficient to reduce the probability of incorrect instantaneous ambiguity fixing by several orders of magnitude.
Finally, the proposed method was tested with three low-cost Multi-GNSS receivers and a magnetic field sensor in a test drive. A tightly coupled Multi-GNSS/INS system served as a reference. The measurement results showed that the calibrated and filtered heading differed by less than 10 • from the reference heading.