Magnetometer Calibration for Small Unmanned Aerial Vehicles Using Cooperative Flight Data

This paper presents a new method to improve the accuracy in the heading angle estimate provided by low-cost magnetometers on board of small Unmanned Aerial Vehicles (UAVs). This task can be achieved by estimating the systematic error produced by the magnetic fields generated by onboard electric equipment. To this aim, calibration data must be collected in flight when, for instance, the level of thrust provided by the electric engines (and, consequently, the associated magnetic disturbance) is the same as the one occurring during nominal flight operations. The UAV whose magnetometers need to be calibrated (chief) must be able to detect and track a cooperative vehicle (deputy) using a visual camera, while flying under nominal GNSS coverage to enable relative positioning. The magnetic biases’ determination problem can be formulated as a system of non-linear equations by exploiting the acquired visual and GNSS data. The calibration can be carried out either off-line, using the data collected in flight (as done in this paper), or directly on board, i.e., in real time. Clearly, in the latter case, the two UAVs should rely on a communication link to exchange navigation data. Performance assessment is carried out by conducting multiple experimental flight tests.


Introduction
The use of unmanned aerial vehicles (UAVs) has been increasing exponentially in recent years, and thus they are having a strong economic impact on the global market [1]. Indeed, they can be used for a large variety of military and civilian applications [2], such as border patrol [3], search and rescue [4], infrastructure monitoring [5], package delivery [6] and precision agriculture [7]. To fully enable the abovementioned mission scenarios, it is widely agreed that strong research efforts must be dedicated to significantly enhance the level of autonomy of UAVs regarding their guidance, navigation and control functions [8], on the one hand, and on the other hand, to the definition of a framework where they can be safely integrated in the civil airspace (see for example the NASA Unmanned Traffic Management program [9] and Corus project by SESAR [10]).
In this respect, this work lies in the framework of research activities, carried out at the Department of Industrial Engineering of the University of Naples "Federico II", aiming at investigating the possibility to exploit cooperative strategies to improve navigation performance of mini/micro UAVs (MAVs) [11,12]. If these vehicles fly under nominal Global Navigation Satellite System (GNSS) coverage, autonomous and safe navigation is enabled by integrating measurements from low-cost GNSS receivers and commercial-grade micro-electro-mechanical systems (MEMS)-based inertial (i.e., gyroscopes In order to address this problem, this paper investigates the possibility to improve the accuracy in the estimate of the magnetic heading provided by magnetometers. The calibration of magnetic sensors is typically based on the assumption that the intrinsic error sources can be characterized using an ellipsoidal error model [22,23]. Consequently, the parameters of the model ellipsoid can be computed by applying a non-linear optimization process, which requires in input a set of measurements taken with different pointing conditions of the sensor's axes in order to properly sample the ellipsoid surface. However, these calibration procedures cannot cope with systematic errors induced by onboard magnetic sources which act during nominal flight operations of the MAV [24,25]. A closed-form analytical solution for the calibration problem of magnetometers on board UAV, conceived to potentially account also for the contribution of the electric motors, can be found in [26]. However, this method has not been tested in flight, and no information is provided by the authors regarding the achievable magnetic heading accuracy. For this reason, an original calibration approach to estimate magnetic biases caused by onboard disturbances is proposed in this paper. This method exploits cooperation between two MAVs, and it requires calibration data to be collected in flight executing a predefined coordinated maneuver. In the following, the vehicle whose magnetometers need to be calibrated is called "chief", while the cooperative vehicle is called "deputy". The chief must be able to detect and track the deputy using a monocular camera (and, consequently, proper vision-based algorithms [27,28]) to estimate the unit vector corresponding to the chief/deputy relative position (i.e., the so-called line-of-sight, LOS) in the camera reference frame. Also, during the calibration flight, both the vehicles must be kept under nominal GNSS coverage to allow estimating their relative position in local (North-East-Down, NED) coordinates. Using the chief/deputy relative position information in camera and NED coordinates, the calibration problem can be mathematically formulated as a system of non-linear equations as described in detail in Section 2. This problem is addressed using a customized implementation of the Levenberg-Marquardt (LM) method [29] (i.e., a state-of-the-art least-squares solver) based on the formulation proposed by Gavin [30]. It is worth outlining that if both the MAVs are equipped with a Sensors 2020, 20, 538 3 of 17 camera, this approach can be used to calibrate the magnetic sensors on board both the two vehicles, before being able to carry out their cooperative or independent mission. Clearly, although this method is conceived to exploit cooperation between two MAVs, the proposed concept is perfectly applicable even if the deputy is another type of autonomous vehicle (e.g., ground or marine) or even if it is a fixed ground station, provided that it can acquire GNSS information while being visually tracked. Finally, it is also important to underline that the calibration approach can be applied either off-line (as done in this paper), or directly on board, i.e., processing in real time the acquired flight data. Clearly, the applicability of this latter option requires that the two UAV are able to exchange information about their navigation status using a communication link. An experimental campaign of flight tests is carried out to assess the applicability and performance of the proposed cooperative calibration method. The flight tests are carried out using both fixed and flying deputies. Preliminary results of this work were initially presented in a conference paper [31].
The rest of the paper is organized as follows: Section 2 formulates the magnetic bias calibration problem, and it describes in detail the proposed cooperative strategy and the associated LM-based algorithmic solution. Section 3 presents the setup used for data collection and the adopted experimental strategy, and it discusses the related results. Finally, Section 4 provides a conclusion and indications about future works.

Magnetic Heading Definition
The calibration problem addressed in this work relies on the assumption that magnetic sensor's measurements are not affected by disturbances from the environment where the MAV is operating. Clearly, this condition is not met if the MAV must fly within areas characterized by large external magnetic sources (e.g., close to power plants, power-lines or other similar facilities). Under this assumption, the magnetic field measured by a three-axes magnetometer (H) can be expressed using Equation (1): where H x , H y and H z are its components in BRF, H E is the Earth's magnetic field, and ∆H is the magnetic bias. ∆H is caused by residual intrinsic error sources and (mostly) by onboard magnetic disturbances, thus being constant in BRF. If the attitude of an MAV (namely, the orientation of its BRF with respect to NED) is parametrized by a 321 sequence of Euler Angles, i.e., heading (ψ), pitch (θ) and roll (ϕ), the Body-Stabilized Reference Frame (BSRF) can be defined as the reference frame obtained by projecting the BRF axes in the North-East plane thanks to the pitch and roll angle estimates provided by the chief onboard autopilot. Thus, the projection of H in the BSRF (H s ) can be computed using Equation (2): where M θ and M ϕ (i.e., the elemental rotation matrixes corresponding to θ and ϕ, respectively) can be written as follows: Sensors 2020, 20, 538 4 of 17 Clearly, the in-plane components of H s (namely, Hx s and Hy s ) can be used to compute the magnetic heading (ψ M ), i.e., the heading angle estimate obtained using the measurements of a three-axes magnetometer, as shown by Equation (4) where δ M is the local magnetic declination, i.e., the angle (on the local horizontal plane) between the magnetic North (N M ) and the true (geographic) North (N). Unfortunately, due to the presence of ∆H, such heading estimate differs from the true heading (ψ) as it is referred to an apparent geographic North direction (N a ) rather than to the true one, as depicted by Figure 1. By substituting Equation (1) in Equation (4), ψ M can be expressed using Equation (5): where H E,xs (∆Hx S ) and H E,ys (∆Hy S ) are the in-plane components of H E (∆H) in BSRF. Consequently, the proposed calibration approach aims at increasing the magnetic heading accuracy by computing the in-plane components of the magnetic bias expressed in BSRF (∆H s ). Indeed, this allows obtaining a calibrated magnetic heading (ψ M,c ) as shown by Equation (6): Clearly, the in-plane components of Hs (namely, Hxs and Hys) can be used to compute the magnetic heading (ψM), i.e., the heading angle estimate obtained using the measurements of a threeaxes magnetometer, as shown by Equation (4 where δM is the local magnetic declination, i.e., the angle (on the local horizontal plane) between the magnetic North (NM) and the true (geographic) North (N). Unfortunately, due to the presence of ∆H, such heading estimate differs from the true heading (ψ) as it is referred to an apparent geographic North direction (Na) rather than to the true one, as depicted by Figure 1. By substituting Equation (1) in Equation (4), ψM can be expressed using Equation (5): where HE,xs (ΔHxS) and HE,ys (ΔHyS) are the in-plane components of HE (ΔH) in BSRF. Consequently, the proposed calibration approach aims at increasing the magnetic heading accuracy by computing the in-plane components of the magnetic bias expressed in BSRF (ΔHs). Indeed, this allows obtaining a calibrated magnetic heading (ψM,c) as shown by Equation (6) Figure 1 provides a simplified representation of the problem addressed in this paper.

Cooperative Calibration Strategy
The determination of ∆Hx S and ∆Hy S must rely on data collected in flight, i.e., when all the chief onboard systems operate nominally, thus ensuring that the level of the associated magnetic disturbances is also nominal. To carry out this task, a cooperative calibration strategy is here proposed. Specifically, the two vehicles must execute a coordinated maneuver, during which the visual camera mounted on board the chief is used to detect and track the deputy through a sequence of frames. Two conditions must be strictly met regarding the definition of the flight path for the two vehicles: • They must fly under nominal GNSS coverage.

•
The deputy must stay within the camera Field-of-View (FOV).
Besides these mandatory points, some additional guidelines can be followed regarding the flight path definition to improve accuracy in the result of the calibration process, as better explained in Section 3.
For each frame acquired by the camera on board the chief, the position of the deputy projection on the image plane can be detected. To this aim, if the magnetic calibration is done in real time, an autonomous visual detection and tracking algorithm must be implemented on board. Due to the cooperative nature of the proposed calibration strategy, this image processing task can be carried out exploiting the knowledge of the deputy (e.g., in terms of size, shape and color) as well as its possibility to exchange navigation data with the chief by means of a communication link, as proposed in [27] and [28]. Instead, if the calibration is carried out off-line (i.e., after flight data acquisition), the visual detection task can be realized using a supervised approach (i.e., the position of the deputy projection on the image plane is extracted manually by a human operator). In both the cases, the image position information can be used to compute the unit vector (u CRF ) representing the deputy LOS in the camera reference frame (CRF) using the intrinsic camera calibration parameters [32]. If autonomous visual algorithms are used, the angular accuracy in the LOS estimate can be of the order of the camera instantaneous FOV [28], while a supervised (off-line) approach can even lead to better accuracy due to the possibility to identify the deputy position at sub-pixel level.
Since the two vehicles are both flying under nominal GNSS coverage, the relative position vector and, consequently, the unit vector representing the deputy LOS in NED with respect to the chief (u NED ) can also be determined. In this respect, it is worth outlining that the level of accuracy characterizing the estimates of u NED depends on the adopted relative positioning method. Specifically, the chief-deputy relative position vector can be obtained either computing the difference between the two position fix or applying Differential Global Positioning System (DGPS) [33], or Carrier-Phase DGPS (CDGPS) [34] algorithms. Clearly, in the latter case, GNSS raw data, namely pseudo-range and carrier-phase measurements, must be available. With regards to the guidelines for flight path definition mentioned earlier, the angular accuracy in u NED can be also improved by keeping the deputy as far as possible from the chief. Of course, the maximum allowable distance is determined by the size of the deputy, the camera resolution and the local background (which can be more or less cluttered) in the image plane.
At this point, if the visual and positioning data are properly synchronized (meaning that a GNSS-based relative position measurement is associated to each visual frame), the relationship between u CRF and u NED can be written as follows: where M M is the camera-body extrinsic calibration matrix and M ψθϕ is the UAV attitude matrix (i.e., the full rotation matrix corresponding to the 321 sequence of Euler angles defined in Section 2.1). If the camera is mounted according to a strapdown configuration, M M (representing the orientation of the CRF with respect to the BRF) is constant and it can be computed by carrying out an ad-hoc extrinsic calibration procedure either off-line [35] or on-line [36]. Using the pitch and roll angles available from the chief onboard navigation system, M ψθϕ can be expressed as a function of the unknown calibrated magnetic heading as shown by Equation (8): Consequently, if Equations (6) and (8) are substituted within Equation (7) for k time instants at which synchronized visual and GNSS data have been acquired, the magnetic bias calibration process can be formulated as a system of 3k non-linear equations in the two unknowns, i.e., ∆Hx S and ∆Hy S . In fact, for each of the k time instants, a triplet of equations is obtained by setting to zero the residual vector function (r) as shown by Equation (9): Hence, the estimation of the unknown magnetic bias components is entrusted to a customized implementation of the LM algorithm, which is an iterative technique to solve non-linear least squares problems. In the following, the symbol r is used to indicate directly the [3k × 1] residual vector obtained writing (9) for all the selected images. A block diagram which summarizes the proposed cooperative calibration approach is depicted in Figure 2.
Consequently, if Equations (6) and (8) are substituted within Equation (7) for k time instants at which synchronized visual and GNSS data have been acquired, the magnetic bias calibration process can be formulated as a system of 3k non-linear equations in the two unknowns, i.e., ΔHxS and ΔHyS. In fact, for each of the k time instants, a triplet of equations is obtained by setting to zero the residual vector function (r) as shown by Equation (9): Hence, the estimation of the unknown magnetic bias components is entrusted to a customized implementation of the LM algorithm, which is an iterative technique to solve non-linear least squares problems. In the following, the symbol r is used to indicate directly the [3k × 1] residual vector obtained writing (9) for all the selected images. A block diagram which summarizes the proposed cooperative calibration approach is depicted in Figure 2.

LM-Based Solution
The LM algorithm allows formulating the least-squares minimization of the residual vector function in Equation (9) as a scalar problem by introducing an equivalent scalar cost function, namely the chi-squared error (χ 2 ). This cost function is the sum of the weighted squares of the errors (also called residuals) between the measured data, i.e., the left-hand side in Equation (7) in this case, and the curve-fit function, i.e., the right-hand side in Equation (7) in this case, and it can be defined as follows:

Hx Hy r Wr
where W is a [3k × 3k] diagonal weight matrix, whose elements are the reciprocal of the weights representing the level of confidence associated to each measured vector (uCRF). Here, W is set to be an identity matrix, meaning that the same weight/importance is assigned to each visual estimate of the deputy LOS in CRF. This choice is in line with the fact that the accuracy in the estimate of the position of the deputy projection on the image plane does not vary significantly with factors like relative

LM-Based Solution
The LM algorithm allows formulating the least-squares minimization of the residual vector function in Equation (9) as a scalar problem by introducing an equivalent scalar cost function, namely the chi-squared error (χ 2 ). This cost function is the sum of the weighted squares of the errors (also called residuals) between the measured data, i.e., the left-hand side in Equation (7) in this case, and the curve-fit function, i.e., the right-hand side in Equation (7) in this case, and it can be defined as follows: where W is a [3k × 3k] diagonal weight matrix, whose elements are the reciprocal of the weights representing the level of confidence associated to each measured vector (u CRF ). Here, W is set to be an identity matrix, meaning that the same weight/importance is assigned to each visual estimate of the deputy LOS in CRF. This choice is in line with the fact that the accuracy in the estimate of the position of the deputy projection on the image plane does not vary significantly with factors like relative distance and local background (which can mostly influence the detection probability, if the image processing task is carried out on line, rather than its accuracy [27]). In this work, the minimization of χ 2 is carried out exploiting the numerical formulation of the LM technique proposed by [30]. According to this formulation, an updated estimate of the unknown in-plane components of ∆H in BSRF is obtained at each LM iteration using Equation (11), where j is the iteration counter, and h LM , i.e., the LM-based correction vector, is computed using the Marquardt's relationship [29]: where λ is the LM damping parameter, J is the [3k × 2] Jacobian matrix defined as in Equation (13), and D JTWJ is a [2 × 2] matrix containing only the diagonal elements of J T WJ: Considering that the LM algorithm is a combination of two minimization techniques, namely the gradient descent method and the Gauss-Newton method, the value of the damping factor is a key parameter to determine its convergence. In fact, for large values of λ, the correction vector obtained from Equation (12) is closer to the one provided by the gradient descent method, meaning that χ 2 is reduced in the direction opposite to its gradient. Instead, for small values of λ, the correction vector obtained from Equation (12) is closer to the one provided by the classical Gauss-Newton method.
Once the LM-based correction is computed, the updated solution provided by Equation (11) is not passed at the next step of the iterative process unless a specific condition on the associated variation on χ 2 is met. The adopted confirmation criterion is given by the following relation [30]: where ε o is a positive threshold set by the user. For instance, this condition allows discarding an update which would cause an increase in χ 2 . Based on the result of Equation (14), the value of the damping factor is changed at each iteration to accelerate convergence. Specifically, if the condition in Equation (14) is met (not met), λ is decreased (increased) using the scaling factor λ DN (λ UN ). Clearly, large scaling factors can provide convergence in fewer iterations, but the process may encounter numerical instabilities. With regards to the initialization of the iterative procedure, the initial values of ∆Hx S and ∆Hy S are set to zero, though a non-zero initial guess may be used if available, e.g., obtained by an on-ground calibration. The initial value of λ (λ 0 ) is also freely selected by the user.
Finally, the iterative process is ended if at least one out of the three criteria shown in Equation (15) is met, unless the iteration counter reaches a maximum value (e.g., 100): It is also worth highlighting that although all the (dimensionless) user parameters, i.e., ε 0 , ε 1 , ε 2 , ε 3 , λ 0 , λ UN and λ DN , defined earlier in this section are common to any LM implementation, their selection must be done in view of the specific application.

Experimental Flight Test Campaign
An experimental campaign of flight tests has been carried out to validate and assess the performance of the proposed cooperative strategy for the calibration of magnetic sensors installed on board MAVs. Specifically, the main performance parameter is the accuracy level charactering the estimation of the calibrated magnetic heading obtained applying Equation (6). To evaluate this accuracy, a reference estimate of the heading angle must be derived to be used as a ground truth (as described in detail in the following). This section is divided in two parts corresponding to different experimental tests carried out using, as deputy, a fixed ground station and a flying vehicle, respectively.

Experimental Setup
For the first experiment, the chief is a customized version of the Pelican TM quadrotor from Ascending Technologies TM (Krailling, GermanY; http://www.asctec.de/) while the role of the deputy is played by a fixed ground station, i.e., a laptop connected via serial RS-232 and a serial-USB adapter to a GNSS antenna (AV59) and receiver (BD960) both produced by Trimble TM (Sunnyvale, CA, USA; https://www.trimble.com/). The chief UAV and the deputy ground station are depicted in Figure 3.
selection must be done in view of the specific application.

Experimental Flight Test Campaign
An experimental campaign of flight tests has been carried out to validate and assess the performance of the proposed cooperative strategy for the calibration of magnetic sensors installed on board MAVs. Specifically, the main performance parameter is the accuracy level charactering the estimation of the calibrated magnetic heading obtained applying Equation (6). To evaluate this accuracy, a reference estimate of the heading angle must be derived to be used as a ground truth (as described in detail in the following). This section is divided in two parts corresponding to different experimental tests carried out using, as deputy, a fixed ground station and a flying vehicle, respectively.

Experimental Setup
For the first experiment, the chief is a customized version of the Pelican TM quadrotor from Ascending Technologies TM (Krailling, GermanY; http://www.asctec.de/) while the role of the deputy is played by a fixed ground station, i.e., a laptop connected via serial RS-232 and a serial-USB adapter to a GNSS antenna (AV59) and receiver (BD960) both produced by Trimble TM (Sunnyvale, CA, USA; https://www.trimble.com/). The chief UAV and the deputy ground station are depicted in Figure 3. The chief uses a 752 × 480 miniaturized CMOS forward looking camera, i.e., the mvBlueFOX-MLC-200wc produced by Matrix Vision TM (Oppenweiler, Germany; https://www.matrixvision.com/), to image the deputy. Since the Pelican autopilot does not give access to the GNSS raw data (which are required to enable DGPS or CDGPS processing), an auxiliary antenna and GNSS receiver (LEA 6T) produced by uBlox TM (Thalwil, Switzerland; https://www.u-blox.com/en) are installed on the chief. Inertial and magnetic sensor measurements are collected at a frequency of 40 Hz, while images and GNSS data are gathered simultaneously at a frequency of 1 Hz. To ensure the synchronization of GNSS, visual and inertial/magnetic data, they must be referred to the same time scale. This task is achieved by using an acquisition software [11] developed to save all the sensor data with an accurate time-tag based on the CPU clock. This time-tag is associated with GNSS measurements, which include the GPS time, with very small latency. Moreover, a precise visual/GNSS synchronization is attained since the image acquisition is triggered by the reception of the first GPS package. The deputy ground station also collects GNSS raw data at 1 Hz.
Despite a single deputy suffices for the application of the proposed calibration approach, an additional ground station is positioned on the test field to collect GNSS raw data. Indeed, if both the The chief uses a 752 × 480 miniaturized CMOS forward looking camera, i.e., the mvBlueFOX-MLC-200wc produced by Matrix Vision TM (Oppenweiler, Germany; https://www.matrix-vision.com/), to image the deputy. Since the Pelican autopilot does not give access to the GNSS raw data (which are required to enable DGPS or CDGPS processing), an auxiliary antenna and GNSS receiver (LEA 6T) produced by uBlox TM (Thalwil, Switzerland; https://www.u-blox.com/en) are installed on the chief. Inertial and magnetic sensor measurements are collected at a frequency of 40 Hz, while images and GNSS data are gathered simultaneously at a frequency of 1 Hz. To ensure the synchronization of GNSS, visual and inertial/magnetic data, they must be referred to the same time scale. This task is achieved by using an acquisition software [11] developed to save all the sensor data with an accurate time-tag based on the CPU clock. This time-tag is associated with GNSS measurements, which include the GPS time, with very small latency. Moreover, a precise visual/GNSS synchronization is attained since the image acquisition is triggered by the reception of the first GPS package. The deputy ground station also collects GNSS raw data at 1 Hz.
Despite a single deputy suffices for the application of the proposed calibration approach, an additional ground station is positioned on the test field to collect GNSS raw data. Indeed, if both the ground stations are contained in the field of view of the camera on board the chief during the entire flight test, the DGPS-vision algorithm described in [11] can be adopted to obtain the desired reference solution for the heading angle of the chief. This quantity, indicated as ψ DGPS-VISION in the following, has an average accuracy of the order of 0.5 • (which is far better than the heading accuracy provided by the Pelican onboard autopilot) [11].
In this work, the proposed magnetic calibration strategy is carried out using off-line the data collected in flight. Of course, to enable an autonomous real-time implementation, a Vehicle-2-Vehicle (V2V) communication link is required to exchange the navigation data. In this respect, it is worth outlining that the amount of information to be exchanged (i.e., the GNSS data) is limited to a few Kbits per second, thus being compatible with the specifications of standard V2V communication systems [37].

Results
The magnetic calibration strategy is applied over data collected during a single flight test, lasting around 8 minutes (although this relatively long time duration is not strictly required for the applicability of the method). The flight path travelled by the chief on the testing site is depicted in Figure 4, where the positions of the two ground stations are also indicated. Among them, the farthest one with respect to the chief (blue dot in Figure 4) is selected as deputy for the cooperative calibration process to limit the uncertainty in the GNSS-based estimation of u NED . The other ground station (green dot in Figure 4) is indicated as auxiliary deputy as it is used to allow the estimation of the reference heading solution (ψ DGPS-VISION ) according to the DGPS-vision algorithm in [11].
flight test, the DGPS-vision algorithm described in [11] can be adopted to obtain the desired reference solution for the heading angle of the chief. This quantity, indicated as ψDGPS-VISION in the following, has an average accuracy of the order of 0.5° (which is far better than the heading accuracy provided by the Pelican onboard autopilot) [11].
In this work, the proposed magnetic calibration strategy is carried out using off-line the data collected in flight. Of course, to enable an autonomous real-time implementation, a Vehicle-2-Vehicle (V2V) communication link is required to exchange the navigation data. In this respect, it is worth outlining that the amount of information to be exchanged (i.e., the GNSS data) is limited to a few Kbits per second, thus being compatible with the specifications of standard V2V communication systems [37].

Results
The magnetic calibration strategy is applied over data collected during a single flight test, lasting around 8 minutes (although this relatively long time duration is not strictly required for the applicability of the method). The flight path travelled by the chief on the testing site is depicted in Figure 4, where the positions of the two ground stations are also indicated. Among them, the farthest one with respect to the chief (blue dot in Figure 4) is selected as deputy for the cooperative calibration process to limit the uncertainty in the GNSS-based estimation of uNED. The other ground station (green dot in Figure 4) is indicated as auxiliary deputy as it is used to allow the estimation of the reference heading solution (ψDGPS-VISION) according to the DGPS-vision algorithm in [11]. The flight path is assigned to provide a significant variation of the heading of the chief. This aspect is important to ensure an adequate sampling of the residual vector function for the minimization process. The time variation of the non-calibrated magnetic heading, obtained applying (4), is depicted in Figure 5. For this computation, the local magnetic declination is obtained using an on-line free software [38], which requires in input the location of the test site (latitude, longitude and altitude) and the time of the experiment. In this case, δM is 3.07°. The flight path is assigned to provide a significant variation of the heading of the chief. This aspect is important to ensure an adequate sampling of the residual vector function for the minimization process. The time variation of the non-calibrated magnetic heading, obtained applying (4), is depicted in Figure 5. For this computation, the local magnetic declination is obtained using an on-line free software [38], which requires in input the location of the test site (latitude, longitude and altitude) and the time of the experiment. In this case, δ M is 3.07 • .
To highlight the effect caused by the presence of onboard magnetic disturbances on the magnetic heading estimation, the heading error (ψ ERR ), measured as the difference between ψ DGPS-VISION and ψ M , is depicted in Figure 6. Specifically, ψ ERR is expressed as a function of time and of the reference ψ DGPS-VISION estimate in Figure 6a,b, respectively.
The starting time in Figure 6a is 20 s. In fact, the time interval during the power on of the electric engines and the takeoff can be neglected as it does not correspond to a nominal flight condition (and, consequently, nominal value of the magnetic bias caused by onboard sources). It is worth outlining that the presence of gaps in the time variation of ψ ERR in Figure 6a is due to temporary absence of the reference heading solution (ψ DGPS-VISION ), mainly occurring if the ground stations are not both enclosed in the camera FOV. By looking at both Figure 6a,b, it is clear how the time variation of ψ ERR is correlated to the variation of the true heading. This is due to the magnetic bias vector (∆H) being constant in BRF. Specifically, the sign of the magnetic heading estimation error is positive for negative heading, while it becomes negative as the chief is pointed with a positive heading.  Both for the computation of ψ M as shown in Figure 5, and for the calibration procedure, the raw magnetic sensor measurements must be filtered to cope with the associated high level of noise. Here, this is done by applying a moving average operator to these noisy magnetic data. A smoothing time of 2.5 s is set when applying the moving average operator considering the relatively smooth flight dynamic and the low magnetometer bandwidth. The result of this operation for the components of H along the x and y axes of BRF is shown by Figure 7.
At this point, the LM-based calibration algorithm can be run. To this end, a key role is played by the selection of the k time instants at which the visual and GNSS data are used. This is done by considering two main guidelines: • Time instants at which the yaw rate, evaluated considering the gyro measurement along the z axis of BRF, is larger than a threshold (τ ω,z ) should be excluded.

•
Time instants at which the chief-deputy relative distance is lower than a threshold (τ R ) should be excluded.
The first condition is set to minimize the heading variation associated to the selected time instants considering the low bandwidth of the magnetometer. The second condition allows limiting the DGPS uncertainty which characterizes the estimation of u NED . In this work, τ ω,z is set to 0.25 • /s, while τ R is set to 70 m, thus allowing the selection of k = 18 time instants. The heading value at the selected frames varies from −38.5 • to 125.3 • , while the mean chief-deputy relative distance is 105.12 m. With regards to the selection of the user's parameters for the LM implementation (defined in Section 2.3), the adopted setting is reported in Table 1. Table 1. Dimensionless parameters characterizing the LM-based magnetic calibration algorithm. The initialization of the damping factor is chosen to favor the Gauss-Newton minimization process with respect to the gradient descent approach, while the values set for λ UN and λ DN have shown to provide good convergence in the LM optimization process [30]. With regards to the convergence thresholds, ε 2 is the largest (10 −4 ) since it depends on the desired accuracy in the estimated value of ∆Hx S and ∆Hy S . In this respect, a variation of ∆Hx S and ∆Hy S in the order of ε 2 must have a negligible effect on the value of r. Given this parameters' setting, the convergence is reached simultaneously for both the first (1 × 10 −11 ) and second (8 × 10 −5 ) condition in Equation (15) after five iterations, while the value of the third convergence metric, i.e., χ 2 /(k−1), stops at (5 × 10 −4 ). The resulting values of ∆Hx S and ∆Hy S are −267.53 and −4.17, respectively.
Using these components of the magnetic bias vector in BSRF, the calibrated magnetic heading can be computed using (6), and the resulting estimation error is reported in Figure 8  By looking at Figure 8b, it can be stated that the proposed calibration approach allows removing the bias induced in the magnetic heading estimation by onboard sources, as the resulting error is independent of ψ DGPS-VISION (unlike in Figure 6b for the non-calibrated magnetic heading). A statistical comparison between the accuracies characterizing ψ M and ψ M,c is reported in Table 2, considering the absolute value of ψ ERR as a performance metric. The statistics reported in the table above quantitatively show the achievable improvement in the magnetic heading accuracy. Specifically, the absolute value of the bias reduces by about 5 • , while the standards deviation reduces by about 1 • .

Experimental Setup
Additional flight tests are conducted using a flying vehicle as deputy for the magnetic calibration process. The advantages of using another UAV as deputy, instead of a fixed ground station, are mainly related to the possibility to realize more quickly the data collection process (while meeting all the identified guidelines regarding the true heading, the range and the yaw rate), especially if the pre-defined coordinated maneuver is executed autonomously. Moreover, as already mentioned, if each UAV is able to visually detect and track the other, the collected data can be used for the calibration of both their corresponding magnetic sensors.
The chief UAV is the same described in Section 3.1.1. Instead, the role of the deputy is played by a X8+ octocopter produced by 3D Robotics TM (Berkeley, CA, USA; https://3dr.com/) depicted in Figure 9. Just like the chief, the deputy is customized by the installation of an auxiliary GNSS receiver and antenna to enable the acquisition of the GNSS raw data. The GNSS receiver is connected via USB to an embedded CPU, i.e., the Odroid XU4 TM now produced by Hardkernel (Anyang, South Korea; https://www.hardkernel.com/), for data acquisition and storage.
Sensors 2020, 20, x FOR PEER REVIEW 13 of 17 to an embedded CPU, i.e., the Odroid XU4 TM now produced by Hardkernel (Anyang, South Korea; https://www.hardkernel.com/), for data acquisition and storage. Besides the deputy, the same ground stations as described in Section 3.1.1 are also installed on the test field and they are always kept in the field of view of the camera to allow the estimation of the reference heading solution. Finally, the same acquisition software and data synchronization strategy described in Section 3.1.1, have been used for these flight experiments. The only difference lies in the fact that inertial data are collected at a higher data rate (i.e., around 76 Hz).

Results
Three flights are carried out at a different test site, so the local value of δM at the time of the experiments is 3.21°. The goal of these additional tests is twofold. On one side, the collected data can be used to further validate the reliability of the proposed method by applying the values of ΔHxS and ΔHyS computed in Section 3.1.2 to correct the bias affecting the estimated magnetic heading. On the other side, the LM-based calibration algorithm can be applied again to check the consistency of the method.
With regards to the first analysis, the time variation of ψERR during the three flight tests for both the calibrated and non-calibrated magnetic heading is reported in Figure 10. Please note that for all the three cases, the time intervals during which the engines are not operating nominally, e.g., before takeoff and after landing, are not considered. Also, during the second and third flight, a set of 360° turns is commanded to the chief, so in the corresponding time interval the ground stations disappear from the FOV and the reference ψDGPS-VISION cannot be computed.  Besides the deputy, the same ground stations as described in Section 3.1.1 are also installed on the test field and they are always kept in the field of view of the camera to allow the estimation of the reference heading solution. Finally, the same acquisition software and data synchronization strategy described in Section 3.1.1, have been used for these flight experiments. The only difference lies in the fact that inertial data are collected at a higher data rate (i.e., around 76 Hz).

Results
Three flights are carried out at a different test site, so the local value of δ M at the time of the experiments is 3.21 • . The goal of these additional tests is twofold. On one side, the collected data can be used to further validate the reliability of the proposed method by applying the values of ∆Hx S and ∆Hy S computed in Section 3.1.2 to correct the bias affecting the estimated magnetic heading. On the other side, the LM-based calibration algorithm can be applied again to check the consistency of the method.
With regards to the first analysis, the time variation of ψ ERR during the three flight tests for both the calibrated and non-calibrated magnetic heading is reported in Figure 10. Please note that for all the three cases, the time intervals during which the engines are not operating nominally, e.g., before takeoff and after landing, are not considered. Also, during the second and third flight, a set of 360 • turns is commanded to the chief, so in the corresponding time interval the ground stations disappear from the FOV and the reference ψ DGPS-VISION cannot be computed. First, it is interesting to note that, since the true heading is positive during the three flight tests, the bias in ψ M is always negative. This is in line with the dependence of the magnetic heading error on the true heading observed in Figure 6b. Overall, Figure 10 shows that the magnetic bias vector components derived thanks to the calibration process described in Section 3.1.2 can be used to correct the magnetometers measurements relative to these new flight tests, thus significantly reducing the systematic error in ψ M,c with respect to ψ M . To synthetize these results, a statistical comparison between the accuracies characterizing ψ M and ψ M,c is reported in Table 3, considering the absolute value of ψ ERR as the performance metric. By looking at these statistics, first, it is important to highlight that the noise in the non-calibrated magnetic heading (measured by the standard deviation of |ψ ERR |) has the same level for all the flight tests analyzed in this work. Instead, the bias (measured by the mean of |ψ ERR |) is larger in the latter three flight tests (i.e., it varies from 11.5 • to 15.2 • ) than in the one analyzed in Section 3.1.2 (i.e., 7.4 • ). This occurs since the constancy of ∆H in BRF induces a bias in the magnetic heading which is a function of the true heading (as highlighted, for instance, in Figure 6b). In the case of the three flights analyzed in this section, the average pointing of the chief (determined by the main direction of the test field with respect to the geographic North) provides a less favorable condition for the vector composition of H E and ∆H (considering that ∆H has its main component along the x axis of BSRF). However, by applying the correction on the magnetic sensor measurements computed by the proposed calibration approach, both the bias and the noise in the estimated magnetic heading are significantly reduced.
Finally, the proposed calibration approach is applied again, but using the data collected during one of the flight tests performed with a flying deputy, in order to show the consistency of the method. In this respect, the third flight is selected as it is characterized by the larger variability in the true heading (which allows better sampling the residual vector function). As τ ω,z and τ R are again set to 0.25 • /s, and 70 m, respectively, the visual and GNSS data required for the calibration are relative to k = 14 time instants. The heading value at the selected frames varies from 47.4 • to 87.2 • , while the mean chief-deputy relative distance is 107.71 m.
By keeping the same LM parameters listed in Table 1, the resulting values of ∆Hx S and ∆Hy S are −275.53 • and 41.16 • , respectively. A statistical comparison between this solution and the one obtained in Section 3.1.2 is carried out in terms of |ψ ERR | computing the calibrated magnetic heading for the third flight test. The results, reported in Table 4, confirm that the same level of accuracy is kept, despite a slight change in the component of ∆H along the y axis of BSRF (which is still much lower than ∆Hx S ).

Conclusions
An innovative method for the calibration of low-cost magnetometers for small/micro unmanned aerial vehicles was presented. It aims at improving the accuracy in the magnetic heading estimate by correcting the systematic error caused by the presence of onboard magnetic fields during nominal flight conditions (e.g., generated by the engines and, potentially, by other electric devices). This problem was addressed by collecting calibration data in flight and adopting a cooperative approach. In fact, the proposed method requires the unmanned aerial vehicle to be calibrated (chief) to visually detect and track a cooperative deputy, both placed under nominal Global Navigation Satellite System coverage to allow relative positioning. If the deputy is a flying vehicle a coordinate flight maneuver must be executed. Using the collected data, the calibration problem can be formulated as a system of non-linear equations solved applying the Levenberg-Marquardt algorithm. An experimental campaign of flight test was conducted to assess the feasibility of the proposed approach and to quantitatively evaluate the achievable magnetic heading accuracy, using either a fixed ground station or a flying vehicle as deputy. Results showed that onboard magnetic sources, being constant in the body reference frame of the chief, induce a systematic error in the estimated magnetic heading, whose entity and sign are a function of the true heading of the chief. This effect is removed by applying the magnetic bias vector components estimated by the calibration procedure to correct the magnetometers measurement. In fact, the residual error in the calibrated magnetic heading is reduced down to two degrees in all the analyzed tests, and it is independent of the true heading of the chief. Besides the systematic error, the noise in the estimated magnetic heading is also slightly reduced after the calibration.
Future research activities will be dedicated to analytically deriving the sensitivity of the calibration approach with respect to the uncertainty in the calibration data (both visual and relative positioning ones). This can also lead to the definition of a standard procedure for data collection which can be efficiently executed exploiting fully autonomous flight mode of the chief and deputy unmanned aerial vehicles. Finally, additional efforts can be addressed to develop a generalized calibration approach able to account also for magnetic disturbances caused by external sources.