1. 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 and accelerometers) and magnetic sensors, within Extended Kalman Filters [
13,
14] or, more recently, Particle Filters [
15], and using either loosely- or tightly-coupled architectures [
16]. In this respect, since onboard gyros do not have enough sensitivity to measure the Earth rate vector, magnetometers play a key role, despite being typically characterized by low bandwidth and high measurement noise. Indeed, their capability to measure the Earth magnetic field, and, consequently to identify the North direction, can be used to bound the heading error by including a magnetic heading estimate among the measurements used by the onboard navigation system of the MAV to get its navigation state. Unfortunately, the presence of magnetic sources, either external or placed on board the MAV, can significantly affect the measurements of a three-axis magnetometer. Together with intrinsic sensor error sources (e.g., biases, scale factor deviations, misalignment and non-orthogonality of the sensor axes) for which a proper ground calibration is always required, these magnetic sources generate disturbances which lead to additional systematic errors in the measurements of the components of the Earth magnetic field in the Body Reference Frame (BRF) of the MAV. Consequently, the measured magnetic heading will be referred to an apparent magnetic North direction, thus leading to an angular heading error which may be up to 10°. For these reasons, such standalone navigation solutions can provide limited accuracy levels, i.e., 5–10 m and 1°–5° for MAVs position and attitude, respectively [
17]. Although, these accuracy levels enable real-time stabilization and control, they do not allow meeting the requirements imposed by several applications (like 3D mapping [
18]), unless ad-hoc solutions are adopted. For instance, multi-antenna-based GNSS techniques allow reaching sub-degree level in the attitude accuracy [
17], but they pose installation constraints which may not be met considering the allocation resources of many micro/small UAVs [
19]. So, in the most general case, considering performance of low-cost IMUs on board MAV, the uncertainties in the roll and pitch estimates provided by the onboard autopilot are typically bounded at the degree-level [
20] (due to the presence of gravity), while the error in the heading angle estimate can be up to several degrees [
21].
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 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.
2. Magnetometers Calibration Method
2.1. 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
Hx,
Hy and
Hz are its components in BRF,
HE 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 (
Hs) 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:
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 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 (
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.
2.2. Cooperative Calibration Strategy
The determination of ΔHxS and ΔHyS 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:
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 (
uCRF) 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 (
uNED) can also be determined. In this respect, it is worth outlining that the level of accuracy characterizing the estimates of
uNED 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
uNED 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
uCRF and
uNED can be written as follows:
where
is the camera-body extrinsic calibration matrix and
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,
(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,
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 3
k 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 [3
k × 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.
2.3. 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
is a [3
k × 3
k] diagonal weight matrix, whose elements are the reciprocal of the weights representing the level of confidence associated to each measured vector (
uCRF). Here,
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
hLM, i.e., the LM-based correction vector, is computed using the Marquardt’s relationship [
29]:
where
λ is the LM damping parameter,
is the [3
k × 2] Jacobian matrix defined as in Equation (13), and
DJTWJ is a [2 × 2] matrix containing only the diagonal elements of
:
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 ΔHxS and ΔHyS 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.
4. 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.