2.1. Nomenclature
Before detailing the adopted methodology, this chapter defines some conventions and symbols that will be used hereafter in the manuscript. Matrices, vectors, and scalars are defined with capital (A), bold (a), and italic (a) font, respectively. The quantity indicates a vector in the frame i. Camera (CRF), Body (BRF), Body Stabilized (BSRF) and locally levelled North East and Down (NED) reference frames are accounted for in this manuscript. Letters c, b, s, and n are used as superscript symbols to indicate CRF, BRF, BSRF and NED frames, respectively. The k-th component of the vector is indicated as , with h = x, y, z. For the sake of clarity, N, D and E subscripts are used when the NED frame is accounted for. A vector measured by an instrument m is indicated by .
The rotation matrix
defines the transformation from frame
i to frame
l, such as
. Attitude angles (heading
ψ, pitch
θ and roll
φ), are defined as the 321-sequence associated to the rotation matrix from NED to BRF, i.e.,
. Indicating with
Mα the elementary rotation of angle α, the NED to BRF rotation matrix can be written as:
Ia and 0a indicate an identity and a zeros matrix of size a × a, respectively. The error on a quantity a is indicated as , which is the difference between the predicted, indicated with a hat (ˆ), and true quantity. The attitude error expressed in NED frame is composed of horizontal, i.e., ρE and ρN, and vertical, i.e., (ρD), components. It connects the predicted and true NED to BRF rotation matrix with , where [×] is the operator yielding the skew symmetric matrix of the vector in the brackets. The Standard Deviation (STD) of scalar and vectorial quantities associated to the measurement instrument m are reported as the scalar σm and the vector Σm, respectively. The derivative of a a × 1 vector a with respect to a b × 1 vector b is a a × b matrix indicated with . The operator that extracts the diagonal from the matrix A is indicated as Diag (A). It returns a vector v, while diag (v) is the operator returning a diagonal matrix, whose elements are the components of the vector v.
2.2. Methodology
The Earth’s magnetic field (
H) components along the BRF directions can be measured by a three-axes magnetometer. From these measurements, the magnetic heading angle can be easily estimated by exploiting the sensed horizontal components of the Earth’s magnetic field as referred to the BSRF, which is obtained as a projection of BRF in the local horizontal north–east plane through a combined pitch (ϑ) and roll (φ) rotation as shown in Equation (2).
The resulting magnetic heading can be therefore expressed as in Equation (3), where
dm is the local magnetic declination, i.e., the angle between the local magnetic (
Nm) and geographical (
N) north directions.
However, onboard electric devices, such as the electric rotors, introduce a disturbance on the quantity sensed by the magnetometers. This disturbance can be modelled as a bias, ΔH, constant in BRF, which induces an error in the estimation of the Earth’s magnetic field direction. Under the assumption of small roll and pitch rotations, ΔH can be assumed to be constant in BSRF, and the effect of the vertical component of ΔH in BRF () on heading estimation can be considered negligible. This assumption is considered in this work since the main interest lies in improving heading and magnetic declination accuracy, without aiming at full 3d calibration.
A graphic illustration of the studied problem can be found in
Figure 1, which clearly shows how the magnetic field vector estimated by magnetometers and projected in BSRF, i.e.,
Hs, can be expressed as the sum of the true Earth magnetic field vector in BSRF (
), which is aligned with the N
m direction, and the internal magnetic biases vector in the same reference frame (Δ
Hs). The effect of such biases can be equivalently seen as an apparent shift of the magnetic and geographic north directions, respectively, to
Nm,a and
Na, as shown in
Figure 1, thus resulting in the heading angle expressed in Equation (3), which will be referred to as “non-calibrated” (
), to be ill-defined. A more accurate angle, referred to as “calibrated” (
), can then be computed if the in-plane components of Δ
Hs are estimated and removed from its formulation, as expressed in Equation (4).
The calibration strategy described in this paper takes its roots from the previous work illustrated in [
18]. As a matter of fact, the bias components
and
are computed by using the same iterative least-squares minimization (LM) procedure. However, the magnetic declination at the flight location is added to the unknowns of the problem. Hence, the proposed approach can be used to correct both “onboard” bias and external disturbances by means of multi-UAV cooperation.
The iterative procedure is based on the minimization of the residual obtained as the difference between the unit vector representing the chief–deputy Line Of Sight (LOS) in NED as reported by CDGNSS data and projected in CRF (
), and the same quantity as estimated by visual-based techniques (
). This can be mathematically put as in Equation (5).
Visual-based detection and tracking techniques [
20] can be used to determine the deputy position in the images collected by the camera on board the chief. Hence, the LOS in CRF is obtained from the pixel coordinates thanks to the intrinsic camera parameters estimated by offline camera calibration. On the other hand, the CDGNSS-based estimate of the LOS in CRF can be obtained by rotating the unit vector corresponding to the CDGNSS relative position vector (
) from NED to BRF, and then from BRF to CRF.
These rotations are achieved by multiplying the vector with the matrices
and
. The misalignment angles between BRF and CRF (determining
), also referred to as extrinsic camera-IMU rotational parameters, can be computed by the camera offline calibration procedure [
21]. More in general, the lever arm of the GNSS antenna with respect to the camera installation on the chief could be considered. However, if the offset is small, the resulting effect is negligible. The final cost vectorial function
r can now be expressed in terms of the unknowns as written in Equation (6).
This calibration technique fully relies on the availability of GNSS data for both platforms and on the deputy visibility with respect to the chief, thus implying that the two vehicles must fly under nominal GNSS coverage conditions and that the deputy must fall in the Field of View (FOV) of the chief camera. For this reason, as a general condition, chief and deputy should be kept facing each other during the flight, if equipped with strapdown frontal-looking cameras. Moreover, the accuracy of the calibration technique is affected by the inter-UAV distance as discussed in detail in
Section 3.1. As a rule of thumb, accurate calibration requires a minimum distance of a few tens of meters.
An additional condition for the proposed calibration procedure, is the need to ensure synchronization of GNSS data (from both the chief and deputy) and camera images, thus building a correspondence between the two vectors used in Equation (6). In this respect, accurate synchronization is performed by time-tagging both chief and deputy data with GNSS time. Clearly, to ensure observability of the problem’s unknowns, Equation (6) must be written considering a set of k camera frames. This results in a 3
k × 1 residual vector
r which can be minimized by applying the LM algorithm. The proposed approach is summarized by the scheme in
Figure 2.
Rather than computing the vectorial cost function as Equation (6) states, an equivalent scalar chi-squared function (χ
2) is computed, and its minimization is carried-out. This function is expressed in Equation (7), where the
W matrix is the 3
k × 3
k block diagonal weight matrix associated to each residual. Such a matrix is shown in Equation (8), where
W1, …,
Wk are the 3 × 3 weight diagonal matrices associated to each of the
k frames and defined, for the
i-th frame, as
. Here,
and
are the covariance matrices associated to CDGNSS and visual measurements expressed in the same reference frame of the residual, i.e., CRF. Detailed derivation of these quantities is reported in
Section 3.
At each new iteration, a correction of the unknown quantities, i.e.,
, is added to their value as computed at the previous step. This quantity is expressed in Equation (9), where
J is the Jacobian matrix representing the derivatives of the cost function
r with respect to the three unknowns identified by the vector
,
λ is the LM damping parameter representing how close the procedure is getting to the gradient descent or the Gauss–Newton methods, and
is the matrix that only contains the diagonal of the
matrix.
The iterative procedure is either stopped when one of the convergence criteria for the LM procedure, as listed in Equation (10), is met or when the maximum number of iterations is reached.