2.1 Model definition

Despite the general idea being known for some time (Koifman and Bar-Itzhack, Reference Koifman and Bar-Itzhack1999), the practical confirmation of VDM-based navigation is relatively recent (Khaghani and Skaloud, Reference Khaghani and Skaloud2018). Indeed, from the large spectrum of possible approaches to apply VDM towards the task of sensor fusion (Sensedobry, Reference Sensedobry2014), only that which uses the VDM as the main process model and inertial and other sensors outputs as observations seems to produce satisfactory results, as investigated in simulations by Cork (Reference Cork2014) for manned aircraft or by Khaghani and Skaloud (Reference Khaghani and Skaloud2016) for a small fixed-wing drone operation at low speeds and limited areas.

In short, this integration scheme uses autopilot's control command (ailerons, elevators, rudder and propeller speed) as an input to the aerodynamic model that specifies forces and moments acting on the drone. Its movement is predicted by resolving equations of motion. The necessary information on the wind is inferred through information redundancy between the aerodynamic model and inertial measurements.

The consistently available data are used as regular observations within the state-space estimation scheme. As detailed by Khaghani and Skaloud (Reference Khaghani and Skaloud2018), the VDM is employed as a system model and the three angular rates and three specific forces of the IMU are used as measurement updates (Khaghani and Skaloud, Reference Khaghani and Skaloud2018, Equation (35)). Data from other sensors are used when available. This is the case for GNSS position and velocity (Khaghani and Skaloud, Reference Khaghani and Skaloud2018, Equations (35) and (36)). The navigation filter (EKF or UKF) estimates the corrections to the navigation ($\boldsymbol {x_n}$) and auxiliary states and its associated covariance matrix ($\boldsymbol {P}$). Due to the large complexity, the elements of the state transition matrix within EKF are derived by an automated linearisation software for the general $l$-frame navigator (Khaghani and Skaloud, Reference Khaghani and Skaloud2018, Equations (7)–(10)), yet its analytical form is presented in full within the Appendix of Khaghani (Reference Khaghani2018) for a local approximation. The same reference contains the models for all formerly mentioned sensors and actuators. The compound error state-vector has the following components and a minimum length:

(1)$$\boldsymbol{x}(46) = [ \boldsymbol{x}_n (12), \, \boldsymbol{x}_p (21), \, \boldsymbol{x}_a (4), \, \boldsymbol{x}_w (3), \, \boldsymbol{x}_e (6) ]^{\rm T},$$

where the auxiliary states are related to sensor (e.g. IMU) time-correlated errors ($\boldsymbol {x}_e$) containing at least three biases (modelled as random constant), actuators errors ($\boldsymbol {x}_a$) containing at least one constant per actuator, wind ($\boldsymbol {x}_w$) corrections per axis as well as VDM model-coefficient errors ($\boldsymbol {x}_p$). The latter is related to all aerodynamic coefficients (hereafter interchangeably referred to as ‘VDM parameters’) quantifying the forces and moments applied to a particular aircraft as depicted in Figure 1 and later detailed in this section. Such platform-dependant coefficients are either only approximately known or they may slightly evolve in time due to some changes in drone payload within a pre-calibrated system. Specific flight patterns and nominal GNSS reception scenarios can be used for such self (in-air) calibration (Laupré and Skaloud, Reference Laupré and Skaloud2021). The approach used to identify IMU stochastic model parameters follows Guerrier et al. (Reference Guerrier, Molinari and Skaloud2015) so $\boldsymbol {x}_e$ may be extended with other time-correlated parameters when required. Their combined effect is subtracted from IMU-specific forces and angular rates before the update. The VDM governing equations following Ducard (Reference Ducard2009) are given in Figure 1 (Equations (4) and (5)). The thrust force $F_T$ is defined in the body frame along the $x$-axis. Here, $V$ is the norm of airspeed i.e. the velocity at which the drone moves through the air; $D$ is the diameter of the propeller; $n$ is the rotation speed of the propeller and advance ratio $J = {V}/{D \pi n}$; $\bar {q}$ is the dynamic air pressure ($\bar {q} = {\rho V^2}/{2}$, with $\rho$ being air density) and $S$ is the wing surface. The moments are expressed in the body frame through its known relation (angles $\alpha$ and $\beta$) to the wind frame. The variables $b$ and $\bar {c}$ correspond to the wingspan and the mean aerodynamic chord of the drone, respectively. The position of the ailerons, the elevators and the rudder are the variables $\delta _a$, $\delta _e$ and $\delta _r$, respectively. Additionally, $\tilde {\boldsymbol {\omega }}=[ b\omega _{l,x},\, \bar {c}\omega _{l,y},\, b\omega _{l,z} ]^{\rm T}/{2V}$, where $\omega _{l}$ is the angular velocity vector of the body frame with respect to the local level frame, correction of which is part of the estimated navigation errors $\boldsymbol {x}_n(12)$, together with position, velocity and attitude.

Figure 1. Left panel shows the used drone with coloured actuating surfaces and relations between airspeed $\boldsymbol {V}$, wind $\boldsymbol {w}$ and drone $\boldsymbol {v}$ velocities through the angle of attack $\alpha$ and side-slip angle $\beta$; right panel shows the modelled moments $M_{xyz}$, thrust force $F_T$ and aerodynamic forces $F_{xyz}$ acting on the drone as a function of respective model coefficients $C_M$, $C_F$, control surface deviations $\delta$, dynamic air-pressure $\bar q$, air-density $\rho$, rotation speed $n$ of propeller with diameter $D$ as well as the advanced ratio $J$ defined in the text together with physical drone parameters $S$, $b$, $\bar c$ and non-dimensional angular rate $\tilde \omega$

In contrast, the platform-dependant geometric parameters such as wing span $b$, mass $m$ and moments of inertia $\boldsymbol {I}^b$ are excluded from $\boldsymbol {x}_p$ since they appear as scaling factors in the model and can be determined with much lower uncertainty a priori as compared to aerodynamic coefficients. Moreover, they are correlated with the aerodynamic coefficients which are already included in $\boldsymbol {x}_p$ (Khaghani and Skaloud, Reference Khaghani and Skaloud2018).

2.2 Model reduction

As previously mentioned, the corrections $\boldsymbol {x}_p$ to aerodynamic coefficients (VDM parameters) are (possibly highly) correlated through the polynomial form of the dynamic model (Equations (4) and (5)). First, such correlations should be considered during filter initialisation with confidence prior to enabling them to be tuned if required. Second, there is a need to combine sets of highly correlated states to reduce the system to a size that is feasible for micro-controller implementation, as was suggested by Laupré and Skaloud (Reference Laupré and Skaloud2020) but not explored so far. Rather than using very well-known approaches of state reduction via sub-optimal filtering (Gelb, Reference Gelb1974), we suggest lumping some correlated parameters for the estimation but using the full model for their application via the re-projected and previously determined correlations. The goal is thus to maintain as much of the system accuracy as possible, yet optimise the system for real-time implementation by suppressing some states, thereby reducing the computational requirements to update the state vector (the number of multiplications increases with the cube of state-vector size). The candidate pairs of coefficients for a possible reduction are identified by analysing the evolution of the covariance matrix $P$ in time. If the correlation between two states is high and does not vary temporally, the updated values are expected to change in a similar way when refined during the state estimation. Linear regressions can be formulated to express one coefficient as a function of the other as

(9)$$C_j = s_{ij} * C_i + o_{ij}$$

where the coefficient $C_j$ is expressed as a scale $s_{ij}$ function of $C_i$ plus an offset $o_{ij}$. The aerodynamic equations can be modified accordingly, reducing the total VDM parameters in the augmented state $\boldsymbol {x}_p$ by the number of candidates found. The choice of the aerodynamic coefficients to be paired and their linear expressions are given in Section 6.2. The influence of the new aerodynamic model on the VDM navigation performance in nominal flight conditions and under GNSS outage are explored in Section 6.6.

3.1 General

The concept of constraining state-space estimators by computer is not new (Gelb, Reference Gelb1974), and for the case of inertial navigation, well analysed (Grewal et al., Reference Grewal, Weill and Andrews2007). Although the general remedies are known (Grewal and Andrews, Reference Grewal and Andrews1993), the potential numerical weaknesses particular to the implementation of VDM-based navigation should be first identified, before deciding on the necessity as well as a strategy for their mitigation. We first recall two types of problems leading to numerical instabilities: (a) ill-conditioning due to large magnitudinal differences in state values (round-off error) and (b) asymmetry of the state covariance matrix $P$. While the symmetry can be forced using a simple yet redundant operation as $P = (P + P^T)/{2}$, round-off errors are complex to avoid and correct as they occur during mathematical operations with limited precision.

An ill-conditioning indicator proposed by Grewal and Andrews (Reference Grewal and Andrews1993) using the largest $\lambda _{{\rm max}}$ and the smallest $\lambda _{{\rm min}}$ eigenvalue of matrix $S$ in the update Kalman filter equations (Table 2) is ${\rm cond}(S)= \lvert \lambda _{{\rm max}} \rvert / \lvert \lambda _{{\rm min}} \rvert$.

Table 2. Summary of discrete Kalman Filter steps

The following rule of thumb can be used to ensure a well-conditioned matrix (Grewal et al., Reference Grewal, Weill and Andrews2007): ${\rm cond}(S) \ll {1}/({2^{-N}})$ with $N$ being the number of bits used in the mantissa. The software environment with a 64-bit architecture uses 52 bits for its mantissaFootnote ¹ and therefore, the condition number should be below $10^{15}$ to guarantee numerical stability of the system. Additionally, round-off errors on the solution of the system $Ax=b$ can be shown to be bound by the following formula (Verhaegen and Van Dooren, Reference Verhaegen and Van Dooren1986):

(10)$$\|(A^{{-}1} b)- (\widehat {A^{{-}1} b}) \| \leq \epsilon \,\cdot\, {\rm cond}(\hat A) \,\cdot\, \|(\widehat {A^{{-}1} b})\|$$

where $\|\,\cdot\, \|$ is the norm, $\epsilon$ is the machine/computer precision, ${\rm cond}(\,\cdot\, )$ is the condition number of a matrix and $(A^{-1} b)$ is the exact solution of the system while the symbol $\hat {(\,\cdot\, )}$ denotes its numerical approximation. The error in the resolution of a linear system is directly proportional to the condition number of the inverted matrix as shown in Equation (10). Although numerical errors do not necessarily lead to filter divergence, the filter can become momentarily unstable, leading to a slower steady-state convergence (Grewal and Andrews, Reference Grewal and Andrews1993) or incorrect state estimation.

3.2 VDM in global-frame

As described by Khaghani and Skaloud (Reference Khaghani and Skaloud2018), the navigation frame is Earth-referenced to the WGS-84 ellipsoid with local level resolving axes and position expressed as geodetic latitude, longitude and height. For numerical operations, latitude and longitude is expressed in $radians$, and height in $m$. When using these units for their corrections within the filter states $\boldsymbol {x}_n$, the values come close to machine precision. Furthermore, the use of differential GNSS decreases the corresponding variance in position. For example, when kept in radians, the achievable centimetre accuracy in position with Real Time Kinematic (RTK) or Post Processing Kinematic (PPK) results in a standard deviation of approximately $10^{-9}$ [rad], which generates in the corresponding covariance matrix $P$ values of $10^{-18}$ [rad$^2$].

Before recalling the need for a re-scaled implementation, the discrete Kalman Filter steps are summarised in Table 2 to define notation, while more details are provided by Grewal and Andrews (Reference Grewal and Andrews1993) and Titterton et al. (Reference Titterton, Weston and Weston2004).

A first step to reduce the condition number is to identify the extreme (smallest and largest) state variances in the initial covariance matrix $\boldsymbol {P}0$: the position ($10^{-18}$ [rad$^2$]) and the propeller speed ($10^2$ [(rad/s)$^2$]). The initial condition number for $P0$ is very high ($10^{22}$), certainly above the aforementioned limit of $10^{15}$. In a second step, chosen scaling factors $S$ are applied to the identified problematic states to reduce the condition number of the covariance matrix $\boldsymbol {P}0$: $S_1 = 10^8$ for latitude/longitude and $S_2 = 10^{-2}$ for the propeller speed. The scaling factors are chosen to change the initial state variances to approach unity. In this work, only the latitude, longitude and propeller speed-related states are scaled ($S_{\varphi }=S_{\lambda }=10^{8}$, $S_n = 0.01$). Practically, a scaled vector $\boldsymbol {x_s}$ is obtained from the initial state vector $\boldsymbol {x}$ multiplied with the scaling diagonal matrix $\boldsymbol {D_x} = [D_1,\, D_2]^{\rm T}$. This transformation can be seen as an arbitrary change of units for (a few) selected variables. Scaling the problematic states needs to be followed by an adaptation of related matrices as detailed in Table 3. As seen in Equation (10), the propagation of round-off errors while performing the matrix inversion is proportional to the condition number. The inversion occurs in the computation of the Kalman gain $\boldsymbol {K}_k$, specifically when the expression $\boldsymbol {G} =\boldsymbol {H} \boldsymbol {P} \boldsymbol {H}^T + \boldsymbol {R}$ is inverted. The units of $\boldsymbol {G}$ are controlled by the measurement noise $\boldsymbol {R}$ and observation $\boldsymbol {H}$ matrices. Hence, in the third step, these matrices must also be scaled with the noise scaling diagonal matrix $\boldsymbol {S}_z$ (if the noise and the states have the same units, $\boldsymbol {S}_z = {\rm diag}(\boldsymbol {S}_x)$ so they correspond to the new state vector $\boldsymbol {x}_s$. We obtain $\boldsymbol {H}_s = \boldsymbol {S}_z \boldsymbol {H}_u\boldsymbol {S}_x^{-1}$ and $\boldsymbol {R}_s = \boldsymbol {S}_z \boldsymbol {R}_u \boldsymbol {S}_z^T$.

Table 3. Adapted scaled matrices

The potentially problematic matrix is thus $\boldsymbol {G}_s = \boldsymbol {S}_z ( \boldsymbol {H}_u \boldsymbol {P}_u {\boldsymbol {H}_{\boldsymbol {u}}}^T +\boldsymbol {R}_u ) \boldsymbol {{S}_{z}}^T$. By carefully choosing the elements of $\boldsymbol {S}_z$, the magnitude of the elements of $\boldsymbol {G}$ can be adjusted to be more homogeneous which in turn lowers the condition number of the matrix. In the same manner, the observations $z$ related to the scaled states should be adapted to match the corresponding scaled observation matrix to compute a meaningful innovation: $\boldsymbol {z} - \boldsymbol {H}\boldsymbol {x} \implies \boldsymbol {S}_{z}\boldsymbol {z} - \boldsymbol {H}_{s}\boldsymbol {x}_{s}$.

In this work, the measurement noise covariance $\boldsymbol {R}$ is chosen to match the sensor characteristics (GNSS and IMU). These are described in Section 5.2. The entries into system noise covariance $\boldsymbol {Q}$ for the relevant states ($\boldsymbol {x}_n$, $\boldsymbol {x}_a$ $\boldsymbol {x}_w$, $\boldsymbol {x}_e$) are presented in Table A1 situated within the Appendix.

4.1 Mixed approach

We revisit two known implementations of the Kalman filter, the combination of which is expected to address: (i) the general challenge of numerical stability of the estimation; (ii) the oscillations of state variables during the in-air initialisation; (iii) the prospect of maintaining the same filtering methodology during nominal and autonomous operation.

The first method exploits one particular form of factorisation put forward by Thornton and Bierman (Reference Thornton and Bierman1978), also known as $UDU$-factorisation, which separates the state covariance matrix as $\boldsymbol {P} = \boldsymbol {U} \boldsymbol {D} \boldsymbol {U}^T$ with $\boldsymbol {U}$ and $\boldsymbol {D}$ an upper triangular and diagonal matrix, respectively. The Kalman filtering steps as presented in Table 2 are adapted accordingly. This factorisation reduces the dynamic ranges of the variables, leading to more homogeneous scales in the computations, and preserves the symmetry of the covariance matrix $P$.

The second method employs the partial Schmidt–Kalman (Schmidt, Reference Schmidt1966) or so-called ‘considered’ filter. The implementation is particularly advantageous when the estimated states are composed in large part of almost constant values such as sensor biases, or in the case of VDM, the aerodynamic coefficients. The partial Schmidt–Kalman filter developed by Brink (Reference Brink2017) is principally defined as a weighted mean between the updated and the predicted ‘partly considered’ states (here abbreviated as ‘considered’ states) and their covariance with a weight factor $\gamma (t) \in \, ]0-1\rangle$ as

(11)$$\begin{align} \hat{\boldsymbol{X}}_{y,k}^c & = \gamma(t) \hat{\boldsymbol{X}}_{y,k} + (1-\gamma(t))\tilde{\boldsymbol{X}}_{y,k} \end{align}$$

(12)$$\begin{align} \hat{\boldsymbol{P}}_{yy,k}^c & = \gamma(t)^2 \hat{\boldsymbol{P}}_{yy,k} + (1-\gamma(t)^2) \tilde{\boldsymbol{P}}_{yy,k} \end{align}$$

The dependency of $\gamma$ on a relative time $t$ is defined later (Equation (14)). The subscript $y$ denotes the elements where the partial Schmidt update is applied. With this definition, $\gamma =1$ computes an optimal Kalman update, whereas $\gamma =0$ creates a Schmidt update on the chosen states.

A practical issue with the $UDU$ decomposition is that this factorisation is not distributive on additions. In a naive implementation, one can reconstruct the covariance matrix, apply the partial reset and factorise it again. However, this would lead to a large increase in computations and partially defeat the purpose, which is to maintain the covariance in a numerically stable form. Carpenter and D'Souza (Reference Carpenter and D'Souza2018) presents a method to apply to a classic Schmidt–Kalman in the case of $UDU$ factorised filters. This development can be adapted to perform the aforementioned weighted mean. Indeed, by using the symmetry of covariance matrices and denoting $(\boldsymbol {H} \boldsymbol {P} \boldsymbol {H}^T + \boldsymbol {R}) = \boldsymbol {G}$, Equation (12) can be reformulated as

(13)$$\begin{align} \hat{\boldsymbol{P}}_{yy,k}^c & = \gamma(t) ^2 ((\boldsymbol{I}-\boldsymbol{K}\boldsymbol{H})\tilde{\boldsymbol{P}}_{k})_{yy} + (1-\gamma(t)^2) \tilde{\boldsymbol{P}}_{yy,k} \nonumber\\ & = (\gamma(t) ^2 \tilde{\boldsymbol{P}}_{k}-\gamma(t) ^2 \boldsymbol{K}\boldsymbol{H}\tilde{\boldsymbol{P}}_{k})_{yy} + (1-\gamma(t)^2) \tilde{\boldsymbol{P}}_{yy,k} \nonumber\\ & =(\tilde{\boldsymbol{P}}_{k}-\gamma(t) ^2 \boldsymbol{K}\boldsymbol{H}\tilde{\boldsymbol{P}}_{k})_{yy} \nonumber\\ & =\hat{\boldsymbol{P}}_{yy,k} + (1- \gamma(t) ^2) (\boldsymbol{K} \boldsymbol{H} \tilde{\boldsymbol{P}}_{k})_{yy} \nonumber\\ & =\hat{\boldsymbol{P}}_{yy,k} + (1- \gamma(t) ^2)(\tilde{\boldsymbol{P}}_{k} \boldsymbol{H} \boldsymbol{G}^{{-}1} \boldsymbol{H} \tilde{\boldsymbol{P}}_{k})_{yy} \nonumber\\ & = \hat{\boldsymbol{P}}_{yy,k} + (1- \gamma(t) ^2)(\tilde{\boldsymbol{P}}_{k} \boldsymbol{H} \boldsymbol{G}^{{-}1} \boldsymbol{G} \boldsymbol{G}^{{-}1} \boldsymbol{H} \tilde{\boldsymbol{P}}_{k})_{yy} \nonumber\\ & =\hat{\boldsymbol{P}}_{yy,k} + (1- \gamma(t) ^2) (\boldsymbol{K} (\boldsymbol{H} \tilde{\boldsymbol{P}}_{k} \boldsymbol{H}^T + \boldsymbol{R}) \boldsymbol{K}^T)_{yy} \end{align}$$

This leads to a similar expression as found by Carpenter and D'Souza (Reference Carpenter and D'Souza2018), with the only difference in the term ($1-\gamma (t)^2$). Since the Bierman–Thornton update step requires scalar processing of measurements, the dimensions of $\boldsymbol {K}$ and $\boldsymbol {G}$ are $n\times 1$ and $1\times 1$, respectively. In other words, an update of rank one is applied directly to the decomposed matrices without the need to recompose them to perform the partial reset. The factorisation adds $n^2$ additions, $n^2 +3n + 2$ multiplications and $n-1$ divisions for each scalar measurement.

Such mixed implementation will be tested outside the nominal scenario within two critical flight phases, namely during the filter initialisation and GNSS signal outage. Both are described in the following subsections. For better clarity, the combination of the $UDU$ factorisation with the partial-Schmidt filter implementation will be referred to abbreviated as ‘partial-Schmidt’ in the following sections.

4.2 Initialisation

At the beginning of the filtering process, the Kalman filter passes through a transient phase. This is partly caused by: (i) the level of uncertainty normally expressed by a quasi diagonal matrix $\boldsymbol {P}0$, as most of the correlations between the states are unknown; (ii) the fact that most of the initial states are unknown (set to zero) and (iii) (some) have large uncertainties relative to the measurement precision. During this phase, there is a high probability that some states converge to a local minimum and remain either incorrectly estimated and/or correlated to other states until the conditions on their observability improve (e.g. due to new dynamics or additional measurements). This in turn limits the navigation quality. The ‘initialisation’ employing a partial Schmidt–Kalman filter is proposed to remedy this problem. To allow the states to evolve more smoothly, these are gradually ‘transited’ from ‘considered’ states to full states. This corresponds to increasing $\gamma$ from 0 to 1 during an initialisation period $\Delta T_{{\rm ini}}$ as

(14)$$\gamma (t) = \dfrac{1}{\Delta T_{{\rm ini}}} (T-T_{{\rm start}}), \quad T \in \langle T_{{\rm start}}, T_{{\rm start}}+\Delta T_{{\rm ini}} \rangle$$

where $T$ is the absolute time in a mission, $T_{{\rm start}}$ is when the initialisation period starts in this time and $t=(T-T_{{\rm start}})$ is the current (relative) estimation time. Such ‘initialisation’ is implemented for the VDM parameters $\boldsymbol {x}_p$ (that are usually pre-calibrated), the wind $\boldsymbol {x}_w$ and the sensor error states $\boldsymbol {x}_e$. These are referenced as ‘considered’ states later on.

4.3 GNSS outage

During a GNSS outage, the absence of position and velocity observations from a GNSS receiver prevents the direct update of navigation states that are related to other auxiliary states via an observation model, resulting in a constant increase in their uncertainties. However, the IMU measurements that are still present, update all navigation states $\boldsymbol {x}_n$ at high frequency. Thus, the elevated variance in position allows for large changes to be applied to the navigation states through IMU updates, possibly leading to erratic jumps in the trajectory. Such corrections are even more exaggerated if the IMU biases are poorly estimated. To counter this effect, the following states are set as ‘considered’ with $\gamma = 0$: the position $\boldsymbol {x}_n(pos)$, the aerodynamic coefficients $\boldsymbol {x}_p$, the wind $\boldsymbol {x}_w$ and the sensor error states $\boldsymbol {x}_e$.

The benefits of the presented combination of filtering strategies during these sensitive phases of flight (i.e. the initialisation and GNSS outage) are presented in Sections 6.4 and 6.5, respectively.

5.1 Platform

The platform used for the experimental flights is depicted in Figure 1(left). Its wing span is $\sim$$1.6$ m and the take-off weight of $\sim$$2.8$ kg. An open source autopilot (FMUv2) from Dronecode (Reference Dronecode2008) was modified to accept the messages from a dual frequency (L1, L2), dual-constellation (GPS, GLONASS) receiver on board (Topcon B110 for the presented tests) and to receive precise pulse per second (PPS) for system-GPS time synchronisation.

The payload is composed of a camera, which is not used within this work, and two IMUs from Thales (Intersence, 2012), installed on a custom micro-controller board (Kluter, Reference Kluter2013), having the following specifications. Gyros: range $\pm$480 deg/s; performance, random walk $0.25$ deg/sqrt(h); noise density $0.004$ deg/s/sqrt(Hz); bias in run stability 12 deg/h; accelerometers: range $\pm$ 8 g; performance, random walk $0.045$ m/s/sqrt(h); noise density, 70 $\mathrm {\mu }$g/sqrt(Hz); bias in run stability $0.1$ mg; sampling, 250 Hz. This board references raw IMU observations in GPS time, and saves and streams them for further processing. Differential post-processing with a reference station allows positioning accuracy to be within a few cm ($0.03$ cm for horizontal, and $0.05$ cm for vertical), velocity accuracy to be within a few cm/s ($0.04$ cm/s for horizontal and $0.08$ cm/s for vertical) and attitude accuracy to be within $\sim$$0.05$–$0.1$ deg. GNSS accuracy is considered time-invariant for simplicity. A Javad Triumph2A (Javad, 2018) is chosen as the GNSS base station and a TOPcon B110 GNSS receiver (TOPcon, 2018) is on the aircraft. The resulting trajectory can be used either for the purpose of calibration and/or reference. The recorded raw inertial, GNSS and autopilot data are replayed as observations for the VDM-based navigation system, the quality of which is assessed with respect to the reference trajectory.

5.2 Flights

Four flights of different geometrical configurations are employed to test the benefit of the proposed numerical strategies within VDM-based navigation on this small fixed-wing drone. First, a 33-minute-long flight is used for calibrating the aerodynamic coefficients and is referred to as flight CF_i8. Three other flights, referred to as AF_i7 (28 min), AF_i6x (23 min) and AF_i6u (17 min) are application flights used to test the suggested modifications to the VDM-based navigation system. As depicted in Figure 2, the flights are dissimilar in their geometry, combining dynamics and a block pattern for CF_i8, different dynamics and block for AF_i7, a long straight corridor for AF_i6x, and a u-shape corridor for AF_i6u. Two flights (CF_i8, AF_i6x) were released as open-source data (Skaloud et al., Reference Skaloud, Cucci and Joseph Paul2021) and we have previously reported them for calibrating the VDM-parameters with the help of attitude updates derived from photogrammetry (Laupré and Skaloud, Reference Laupré and Skaloud2021). Although possible, such calibration is less practical and therefore not used here.

Figure 2. Experimental flights references (blue): (a) CF_i8; (b) AF_i7; (c) AF_i6x and (d) AF_i6u, beginning of the trajectory (red triangle)

5.3 Methodology

At first, the VDM-based navigation system is executed on flight CF_i8 with the nominal (i.e. uninterrupted) reception of GNSS signals for the purpose of self-calibration of VDM coefficients. These are estimated via an optimal recursive smoother. The obtained $\boldsymbol {x_p}$ and corresponding block-covariance matrix $\boldsymbol {P_p}$ are used as priors for all other application flights while increasing the uncertainty of $1\sigma$ to $2\%$ of its initial values. During the three application flights, the VDM parameters are fine-tuned (whenever GNSS positioning is present) to allow for potential small refinements due to changes in the platform configuration and its environment related to battery position, the re-assembled parts with a slightly different orientation between flights, small modification on the payload or changes in the weather conditions.