Unit dual quaternion‐based pose optimisation for visual runway observations

Abstract

This study addresses the pose estimation problem of an aircraft runway using visual observations in a landing approach scenario. The authors utilised the fact that the geodetic coordinates of most runways are known precisely with highly visible markers. Thus, the runway observations can increase the level of situational awareness during the landing approach, providing additional redundancy of navigation and less reliance on global positioning system. A novel pose optimisation algorithm is proposed utilising unit dual quaternion for the runway corner observations obtained from a monocular camera. The estimated runway pose is further fused with an inertial navigation system in an extended Kalman filter. An open‐source flight simulator is used to collect and process the visual and flight dataset during the landing approach, demonstrating reliable runway pose estimates and the improved inertial navigation solution.

Introduction

Pose estimation is a fundamental problem of determining the position and attitude (or orientation) of one object's coordinate system with respect to another. An extension of this principle to computer vision transforms the problem into estimating the pose between the object and the camera from which it is being observed. There has been significant progress in determining the pose of six degrees-of-freedom platforms, particularly by integrating inertial measurement units and monocular/stereo cameras for navigation and mapping [1–5].

For crewed aircraft, studies have shown that the lack of positional awareness is a major cause of accidents [6]. In particular, estimating the pose of a runway relative to a monocular camera offers benefits to all aircraft which should not be underestimated. Moreover, it can be used to reduce errors in the on-board navigation system and adds a level of redundancy to the system. A considerable amount of research has gone into vision-based landing scenarios [7, 8], including the pose estimation for landings on aircraft-carrier ships [9–11]. Most of the work has focused on vertical-take-off-and-landing platforms and not much on the fixed-wing types, except the several recent works such as in [12]. The reason is due to the small margin of the approaching angle in the most ship-landing scenario and thus the difficulty in obtaining precise navigational information from the cameras relative to the runway.

In the application presented here, we consider an airport runway as the object on which the aircraft is to land. Several characteristics of runways make them highly suitable as navigation aids, especially when an aircraft is on its final approach [13, 14]. First, the precise geodetic coordinates of most runways are well known, thus allowing them to be used as absolute navigation references. Secondly, the dimensions of runways must adhere to strict protocols, as does the white markers placed on them. Thirdly, runways are designed to be highly visible, so in essence, they can be used as a navigation tool even if the aircraft does not intend to land at the runway it is observing. Fourthly, there is no limit to the number of aircraft that can observe a given runway for navigation. The emphasis in this work is to increase the level of situational awareness for an aircraft on its final approach.

We propose a novel optimisation algorithm by parameterising the runway pose as a unit dual quaternion (UDQ). The UDQ has been widely applied for the pose estimation problem such as strap-down inertial navigation system [15], coordination of multiple rigid bodies [16], and pose-graph simultaneous localisation and mapping (SLAM) [17]. The UDQ has the advantages of no singularity, unlike Euler angles, and fast convergence due to the unified optimisation of the non-Euclidian rotation and translational vector. It should be noted that the reliable detection of the runway (for example, detecting the piano-key -like markers on the runway) is important and has been the subject of studies [18, 19], but it is not the focus in this work. The unit vectors describing the four runway corners in the camera frame are used to formulate a cost function, which is derived from the geometric relationships between the observations and runway, transformed into the dual quaternion algebra. The pose is estimated through the Levenberg–Marquardt (LM) optimisation method. The contributions of this work are

• UDQ parameterisation to estimate the runway pose from monocular images;

• camera-in-the-loop demonstration of visual-inertial navigation by utilising an open-source flight simulator.

To the author's knowledge, there is no prior literature investigating the UDQ-based optimisation for the runway pose estimation problem. We also utilise a FlightGear simulator to demonstrate the method, as shown in Fig. 1, illustrating an example of image observations of a runway using the simulator, as well as showing the extracted visual features. The simulator is open source and flexible, and it can be effectively used for the camera-in-the-loop simulation, which we believe is a valuable tool for vision-based research.

[IMAGE OMITTED. SEE PDF]

The remaining part of this paper is outlined as follows. Section 2 provides a brief overview of the UDQ algebra, and Section 3 details the UDQ parametrisation for the runway observations obtained from a monocular camera followed by the LM optimisation. Section 4 shows the integration of the runway observations and inertial navigation system. Experimental results are given in Section 5, providing the flight simulator interface, and analysing the optimisation and filtering performance. Section 6 will conclude with future direction.

UDQ algebra

The UDQ is attractive for the estimation of the pose, as the dual quaternion can be represented as a single vector that is amenable to implementation in a statistical estimation process. In addition to this, the rotational and translational parts can be simultaneously optimised in a cost function setting, which is separable and many cases and expandable to include any number of observations. It also has a closed-form solution to the Taylor series expansion [$\mathit{f}(\mathit{q}+ϵ\mathit{p})=\mathit{f}(\mathit{q})+ϵ\dot{\mathit{f}}(\mathit{p})$ for arbitrary quaternions q and p as ${ϵ}^2={ϵ}^3=\cdots =0$].

The dual quaternion comprises two parts that can describe the rigid-body transformation of one coordinate frame with respect to another and is given by 1 $\stackrel{⌣}{\mathit{q}}=\mathit{q}+ϵ\mathit{p}$ where the $\stackrel{⌣}{\mathit{q}}$ is used to signify the dual quaternion, and both q and p are quaternions representing the rotation and translation, respectively, while $\mathit{q}={\mathit{q}}_0+\mathit{i}{\mathit{q}}_1+\mathit{j}{\mathit{q}}_2+\mathit{k}{\mathit{q}}_3:=[{\mathit{q}}_0,{\mathit{q}}_1,{\mathit{q}}_2,{\mathit{q}}_3]$ [Both of the algebraic and vector forms of the quaternion will be used interchangeably in this work.]. The parameter $ϵ$ is used simply to enforce a multiplicative rule onto dual quaternions, being ${ϵ}^2=0$. The product of two dual quaternions is computed as 2 $\begin{array}{rl}{\stackrel{⌣}{\mathit{q}}}_1{\stackrel{⌣}{\mathit{q}}}_2& =({\mathit{q}}_1+ϵ{\mathit{p}}_1)({\mathit{q}}_2+ϵ{\mathit{p}}_2)\\ & ={\mathit{q}}_1{\mathit{q}}_2+ϵ({\mathit{q}}_1{\mathit{p}}_2+{\mathit{p}}_1{\mathit{q}}_2)\end{array}$

3 $=[{\stackrel{⌣}{\mathit{q}}}_1]{}_{+}{\stackrel{⌣}{\mathit{q}}}_2$ 4 $=[{\stackrel{⌣}{\mathit{q}}}_2]{}_{-}{\stackrel{⌣}{\mathit{q}}}_1,$ where $[\stackrel{⌣}{\mathit{q}}]{}_{+}$ and $[\stackrel{⌣}{\mathit{q}}]{}_{-}$ are the left- and right-product matrices of a dual quaternion, respectively, defined as 5 $\begin{array}{rl}[\stackrel{⌣}{\mathit{q}}]{}_{+}:& =\left[\begin{array}{cc}{[\mathit{q}]}_{+}& 0_{4\times 4}\\ {\left[\mathit{p}\right]}_{+}& {[\mathit{q}]}_{+}\end{array}\right]\\ [\stackrel{⌣}{\mathit{q}}]{}_{-}:& =\left[\begin{array}{cc}{[\mathit{q}]}_{-}& 0_{4\times 4}\\ {\left[\mathit{p}\right]}_{-}& {[\mathit{q}]}_{-}\end{array}\right],\end{array}$ in which, $[\mathit{q}]{}_{+}$ and $[\mathit{q}]{}_{-}$ are the standard left- and right-product matrices of a quaternion defined as 6 $[\mathit{q}]{}_{+}:=\left[\begin{array}{cccc}{\mathit{q}}_0& -{\mathit{q}}_1& -{\mathit{q}}_2& -{\mathit{q}}_3\\ {\mathit{q}}_1& {\mathit{q}}_0& -{\mathit{q}}_3& {\mathit{q}}_2\\ {\mathit{q}}_2& {\mathit{q}}_3& {\mathit{q}}_0& -{\mathit{q}}_1\\ {\mathit{q}}_3& -{\mathit{q}}_2& {\mathit{q}}_1& {\mathit{q}}_0\end{array}\right]$ 7 $[\mathit{q}]{}_{-}:=\left[\begin{array}{cccc}{\mathit{q}}_0& -{\mathit{q}}_1& -{\mathit{q}}_2& -{\mathit{q}}_3\\ {\mathit{q}}_1& {\mathit{q}}_0& {\mathit{q}}_3& -{\mathit{q}}_2\\ {\mathit{q}}_2& -{\mathit{q}}_3& {\mathit{q}}_0& {\mathit{q}}_1\\ {\mathit{q}}_3& {\mathit{q}}_2& -{\mathit{q}}_1& {\mathit{q}}_0\end{array}\right].$ Please note that the matrix representations of quaternions and dual quaternions are for the convenience of the product operations and are thus adopted in this work.

Similar to the unit quaternion, a UDQ satisfies the unit norm 8 ${\stackrel{⌣}{\mathit{q}}}^{\ast }\stackrel{⌣}{\mathit{q}}=1$ where ${\stackrel{⌣}{\mathit{q}}}^{\ast }={\mathit{q}}^{\ast }+ϵ{\mathit{p}}^{\ast }$ is the conjugate dual quaternion of $\stackrel{⌣}{\mathit{q}}$. Substituting to (1), we can obtain two intrinsic constraints of the UDQ 9 $\parallel \mathit{q}\parallel =1,{\mathit{q}}^{\ast }\mathit{p}=0.$ It is well known that a UDQ has specific meanings for the rigid-body transformation, which are nicely presented by the following properties.

$$

Property 1

Let $\mathit{t}=[\mathit{x},\mathit{y},\mathit{z}]{}^{\mathrm{T}}$ and q be the translation vector and quaternion of a rigid body with respect to a global frame, respectively, then the rigid-body transformation can be represented as a following UDQ [17]: 10 $\stackrel{⌣}{\mathit{q}}=\mathit{q}+ϵ\mathit{p},\mathrm{with}\mathit{p}=\frac{1}{2}\mathit{t}\mathit{q},$ where $\mathit{t}=[0,\mathit{t}]{}^{\mathrm{T}}$ is the quaternion form of the vector $\mathit{t}$.

Runway observation using UDQ

The pose measurement of the runway $\mathit{z}$ is parametrised in UDQ form $\stackrel{⌣}{\mathit{q}}$ : 11 $\mathit{z}(\mathit{q}+ϵ\mathit{p})=[\mathit{x},\mathit{y},\mathit{z},{\mathit{q}}_0,{\mathit{q}}_1,{\mathit{q}}_2,{\mathit{q}}_3]{}^{\mathrm{T}},$ which are obtained from the on-board vision and inertial sensors. Each of the four corner locations in the image can be transformed into unit direction vectors using the known optical properties of the sensor. This relationship is illustrated in Fig. 2. The objective is to derive the equations that relate the aircraft rotation and translation relative to the runway. The runway state relative to the navigation frame can then be obtained through its relationship with the aircraft. However, the relative rotation and translation of the aircraft relative to the runway must first be ascertained from the corner observations. This is achieved using dual quaternions, and the output of this transformation process is used as the runway observation. The unit directional corner position vectors provide a unique solution to the aircraft-runway attitude. By assuming a pinhole-camera model, each vector is obtained from the $(\mathit{x},\mathit{y})$ coordinates of its respective corner on the image plane and, the focal length of the camera (in pixels). The four unit vectors in the camera coordinate system, denoted $\{{\hat{\mathit{p}}}_{\mathit{i}}^{\mathrm{c}}\in {\mathbb{R}}^3\}{}_{\mathit{i}=\{1,\dots ,4\}}$, are precipitated through a normalisation process.

[IMAGE OMITTED. SEE PDF]

The problem is formulated by choosing any two unit vectors, $\{{\mathit{p}}_{\mathit{i}}^{\mathrm{c}},{\mathit{p}}_{\mathit{j}}^{\mathrm{c}}\}{}_{\mathit{i}\ne \mathit{j}}$, the plane defined by them, and the unit normal vector to this plane, denoted as ${\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}$. It is obtained by taking the cross product of $\{{\mathit{p}}_{\mathit{i}}^{\mathrm{c}},{\mathit{p}}_{\mathit{j}}^{\mathrm{c}}\}$ and normalising the result. Also, required is the unit vector, ${\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}$, that describes the direction of the vector created by subtracting one position vector from the other, relative to the runway frame. From the definitions, it can be observed that ${\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}\in \mathrm{span}\{{\mathit{p}}_{\mathit{i}}^{\mathit{n}},{\mathit{p}}_{\mathit{j}}^{\mathit{n}}\}$. Then, it is necessarily the case that ${\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}\cdot {\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}=0$. This provides the constraint 12 ${\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}\cdot \left({\mathit{R}}_{\mathit{n}}^{\mathrm{c}}{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}\right)=0,$ where the direction vector has been given with respect to the runway navigation frame, and the rotation matrix transforms a vector in the navigation frame to the camera frame: $({\mathit{R}}_{\mathit{n}}^{\mathrm{c}}){}^{\mathrm{T}}(\mathit{q})=\left[\begin{array}{ccc}1-2({\mathit{q}}_2^2+{\mathit{q}}_3^2)& 2({\mathit{q}}_1{\mathit{q}}_2+{\mathit{q}}_0{\mathit{q}}_3)& 2({\mathit{q}}_1{\mathit{q}}_3-{\mathit{q}}_0{\mathit{q}}_2)\\ 2({\mathit{q}}_1{\mathit{q}}_2-{\mathit{q}}_0{\mathit{q}}_3)& 1-2({\mathit{q}}_1^2+{\mathit{q}}_3^2)& 2({\mathit{q}}_2{\mathit{q}}_3+{\mathit{q}}_0{\mathit{q}}_1)\\ 2({\mathit{q}}_1{\mathit{q}}_3+{\mathit{q}}_0{\mathit{q}}_2)& 2({\mathit{q}}_2{\mathit{q}}_3-{\mathit{q}}_0{\mathit{q}}_1)& 1-2({\mathit{q}}_1^2+{\mathit{q}}_2^2)\end{array}\right].$ Our rectangular runway as defined by four points facilitates ${\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}$ expression by utilising each pair of points to define the direction of their connecting edge relative to the runway frame. Thus, the six possible combinations reduce to four vectors for ${\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}$ 13 ${\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}=\left\{\begin{array}{ccc}(1& 0& {0)}^{\mathrm{T}}\\ (0& 1& {0)}^{\mathrm{T}}\\ (({\mathit{x}}_1^{\mathit{n}}-{\mathit{x}}_4^{\mathit{n}})& {\mathit{y}}_1^{\mathit{n}}& {0)}^{\mathrm{T}}\\ (({\mathit{x}}_2^{\mathit{n}}-{\mathit{x}}_3^{\mathit{n}})& {\mathit{y}}_2^{\mathit{n}}& {0)}^{\mathrm{T}}\end{array}\right.,$ where the first and second entries represent the edges (1, 2) (or (3, 4)) and (1, 3) (or (2, 4)), respectively, with the third and fourth acquired from the two diagonal combinations of corners (1, 4) and (2, 3). The second equation is obtained through a transformation and requires the dimensions of the runway to be known. The position of the j th corner, relative to the centre of the runway frame ${\mathit{t}}^{\mathrm{c}}$ is ${\mathit{p}}_{\mathit{j}}^{\mathit{n}}$. This can be expressed relative to the camera frame as 14 ${\mathit{p}}_{\mathit{j}}^{\mathrm{c}}={\mathit{R}}_{\mathit{n}}^{\mathrm{c}}{\mathit{p}}_{\mathit{j}}^{\mathit{n}}+{\mathit{t}}^{\mathrm{c}},$ which must also be perpendicular to ${\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}$, producing a second constraint 15 ${\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}\cdot ({\mathit{R}}_{\mathit{n}}^{\mathrm{c}}{\mathit{p}}_{\mathit{j}}^{\mathit{n}}+{\mathit{t}}^{\mathrm{c}})=0.$ Expanding ${\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}$ and substituting into (12) and (15) produces the two complete homogeneous equations 16 $\left(\frac{{\mathit{p}}_{\mathit{i}}^{\mathrm{c}}\times {\mathit{p}}_{\mathit{j}}^{\mathrm{c}}}{\parallel {\mathit{p}}_{\mathit{i}}^{\mathrm{c}}\times {\mathit{p}}_{\mathit{j}}^{\mathrm{c}}\parallel }\right)\cdot {\mathit{R}}_{\mathit{n}}^{\mathrm{c}}{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}=0$ 17 $\left(\frac{{\mathit{p}}_{\mathit{i}}^{\mathrm{c}}\times {\mathit{p}}_{\mathit{j}}^{\mathrm{c}}}{\parallel {\mathit{p}}_{\mathit{i}}^{\mathrm{c}}\times {\mathit{p}}_{\mathit{j}}^{\mathrm{c}}\parallel }\right)\cdot ({\mathit{R}}_{\mathit{n}}^{\mathrm{c}}{\mathit{p}}_{\mathit{j}}^{\mathit{n}}+{\mathit{t}}^{\mathrm{c}})=0.$ The solving of which produces the position ${\mathit{t}}^{\mathrm{c}}$ and attitude ${\mathit{R}}_{\mathit{n}}^{\mathrm{c}}$ of the camera relative to the runway frame.

In this application, the solutions are obtained by reformulating (16) and (17) using a UDQ, $\stackrel{⌣}{\mathit{q}}=\mathit{q}+\epsilon \mathit{p}$, as the error functions of the LM optimisation. The reformulation process necessitates the transformation of unit vectors and rotation matrices into their quaternion equivalents, in which a vector $\mathit{v}$ changes to a quaternion equivalent form $\mathit{v}=[0,{\mathit{v}}^{\mathrm{T}}]{}^{\mathrm{T}}$, and the rotational operation of a vector ${\mathit{v}}^{\mathit{n}}$ to ${\mathit{v}}^{\mathrm{c}}$ becomes ${\mathit{v}}^{\mathrm{c}}={\mathit{q}}^{\ast }{\mathit{v}}^{\mathit{n}}\mathit{q}$. By utilising the matrix form of the quaternions, the alternative representations of (16) and (17) become 18 $\begin{array}{rl}{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}\cdot \left({\mathit{R}}_{\mathit{n}}^{\mathrm{c}}{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}\right)& =({\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}){}^{\mathrm{T}}\left([\mathit{q}]{}_{+}^{\mathrm{T}}{[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]}_{+}\mathit{q}\right)\\ & ={\mathit{q}}^{\mathrm{T}}\left([{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}{[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]}_{+}\right)\mathit{q}\end{array}$ 19 $\begin{array}{rl}{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}\cdot ({\mathit{R}}_{\mathit{n}}^{\mathrm{c}}{\mathit{p}}_{\mathit{j}}^{\mathit{n}}+{\mathit{t}}_{\mathit{r}})& ={\mathit{q}}^{\mathrm{T}}\left([{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}{[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]}_{+}\right)\mathit{q}\\ & +{\mathit{q}}^{\mathrm{T}}\left(2[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{+}^{\mathrm{T}}\right)\mathit{p},\end{array}$ where p is the dual part of the dual quaternion $\stackrel{⌣}{\mathit{q}}$, given by $\mathit{p}=(1/2)\mathit{t}\mathit{q}$.

The inclusion of the unit-magnitude constraints of the dual quaternion ${\mathit{q}}^{\mathrm{T}}\mathit{q}=1$ and ${\mathit{q}}^{\mathrm{T}}\mathit{p}=0$ (in a vector form) yields the cost functional 20 $\mathit{J}(\stackrel{⌣}{\mathit{q}})=\sum _{\mathit{i}=1}^6\left(\begin{array}{c}{\left({\mathit{q}}^{\mathrm{T}}[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}{[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]}_{+}\mathit{q}\right)}^2\\ +{\left({\mathit{q}}^{\mathrm{T}}[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]\mathit{q}+{\mathit{q}}^{\mathrm{T}}2[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{+}^{\mathrm{T}}\mathit{p}\right)}^2\\ +{\lambda }_1{\left({\mathit{q}}^{\mathrm{T}}\mathit{q}-1\right)}^2+{\lambda }_2{\left({\mathit{q}}^{\mathrm{T}}\mathit{p}\right)}^2\end{array}\right),$ where ${\lambda }_1$ and ${\lambda }_2$ are the Lagrangian multipliers, and the objective is to select values of q and p that minimise the error output, which in this case, is zero. The key property of this cost functional is that the error function is separable, allowing q to be solved independently of p. Its alternate form is given by 21 ${\mathit{J}}_1(\mathit{q})=\sum _{\mathit{i}=1}^6{\left({\mathit{q}}^{\mathrm{T}}[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}{[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]}_{+}\mathit{q}\right)}^2+{\lambda }_1{\left({\mathit{q}}^{\mathrm{T}}\mathit{q}-1\right)}^2$ 22 ${\mathit{J}}_2(\mathit{q},\mathit{p})=\sum _{\mathit{i}=1}^6\left(\begin{array}{c}{\left({\mathit{q}}^{\mathrm{T}}[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]\mathit{q}+{\mathit{q}}^{\mathrm{T}}2[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{+}^{\mathrm{T}}\mathit{p}\right)}^2\\ +{\lambda }_2{\left({\mathit{q}}^{\mathrm{T}}\mathit{p}\right)}^2\end{array}\right),$ where the real-valued quaternion q is the first to be obtained through the LM algorithm. Using the Jacobian $\mathrm{\nabla }{\mathit{J}}_1(\mathit{q})$, the LM algorithm searches in the direction given by the estimated solution $\hat{\mathit{q}}(\mathit{k})$ to the equation 23 $[\mathrm{\nabla }{\mathit{J}}_1^{\mathrm{T}}\mathrm{\nabla }{\mathit{J}}_1+\lambda \mathit{I}]\delta \hat{\mathit{q}}(\mathit{k})=-(\mathrm{\nabla }{\mathit{J}}_1^{\mathrm{T}}){\mathit{J}}_1(\hat{\mathit{q}}(\mathit{k}))$ where $\delta \mathit{q}(\mathit{k})$ is the estimated quaternion error at k -iteration, and $\lambda$ is a non-negative scalar and $\mathit{I}$ is the identity matrix, while the Jacobian matrix $\mathrm{\nabla }{\mathit{J}}_1(\hat{\mathit{q}}(\mathit{k}))$ evaluates to 24 $\mathrm{\nabla }{\mathit{J}}_1(\hat{\mathit{q}}(\mathit{k}))=\sum _{\mathit{i}=1}^6\left(\begin{array}{c}2\left({\hat{\mathit{q}}}^{\mathrm{T}}[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}{[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]}_{+}\hat{\mathit{q}}\right)\\ \times \left({\hat{\mathit{q}}}^{\mathrm{T}}[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathrm{c}}]{}_{-}^{\mathrm{T}}{[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]}_{+}+{\hat{\mathit{q}}}^{\mathrm{T}}[{\mathit{n}}_{\mathit{i}\mathit{j}}^{\mathit{c}}]{}_{-}^{\mathrm{T}}[{\mathit{d}}_{\mathit{i}\mathit{j}}^{\mathit{n}}]{}_{+}^{\mathrm{T}}\right)\\ +4{\lambda }_1\left({\hat{\mathit{q}}}^{\mathrm{T}}\hat{\mathit{q}}-1\right){\hat{\mathit{q}}}^{\mathrm{T}}\end{array}\right).$ During the operation, the estimate $\delta \hat{\mathit{q}}(\mathit{k})$ is used to update $\hat{\mathit{q}}(\mathit{k})$, and the estimate reduces with successive iterations. Convergence is achieved when the magnitude of the cost functional is less than a preset threshold (${10}^{-4}$ is used in this work). Convergence was consistently attained with <18 iterations; the time taken for which was <0.1 s, rendering this method appropriately for real-time implementation. For the generic case, the LM algorithm is guaranteed to converge but may not necessarily converge to the global minimum of the function. In this particular application, however, they are one and the same.

The estimate of the dual part of the quaternion $\mathit{p}(\mathit{k})$ can be obtained by substituting the final attitude estimate $\mathit{q}(\mathit{k})$ into the Jacobian of (20) or (21). Alternatively, it can be obtained through its relationship with $\mathit{q}(\mathit{k})$ 25 $\mathit{p}(\mathit{k})=-({\mathit{A}}_1+{\mathit{A}}_2^{\mathrm{T}}){}^{-1}{\mathit{A}}_2\mathit{q}(\mathit{k}),$ with 26 ${\mathit{A}}_1=\left(\sum _{\mathit{i}}^{\mathit{k}}{\beta }_{\mathit{i}}\right)\mathit{I}$ 27 ${\mathit{A}}_2=2\sum _{\mathit{i}}^{\mathit{k}}{\beta }_{\mathit{i}}([{\mathit{p}}_{\mathit{i}}^{\mathit{n}}]{}_{+}-[{\mathit{p}}_{\mathit{i}}^{\mathit{n}}]{}_{-})\mathit{q}(\mathit{k}).$ where ${\beta }_{\mathit{i}}$ represents the least-square coefficients of the LM minimisation algorithm. The UDQ $\stackrel{⌣}{\mathit{q}}=\mathit{q}+ϵ\mathit{p}$ of the runway is thus obtained as the output of the LM optimisation. If necessary, $\mathit{t}=[\mathit{x},\mathit{y},\mathit{z}]{}^{\mathrm{T}}$ can be computed from the UDQ by utilising (10).

Runway-aided inertial navigation

To demonstrate the effectiveness of the runway pose observation, a runway-aided inertial navigation system is designed. The nominal inertial navigation model is a simplified one utilising a local-fixed, local-tangent navigation frame which is suitable for low-quality inertial sensor applications. An integration extended Kalman filter consists of an inertial state vector $\mathit{x}=[({\mathit{t}}^{\mathit{n}}){}^{\mathrm{T}},({\mathit{v}}^{\mathit{n}}){}^{\mathrm{T}},({\mathit{q}}^{\mathit{n}}){}^{\mathrm{T}}]{}^{\mathrm{T}}$ with a state kinematic model: ${\dot{\mathit{t}}}^{\mathit{n}}={\mathit{v}}^{\mathit{n}}+{\mathit{w}}_1$ ${\dot{\mathit{v}}}^{\mathit{n}}={\mathit{R}}_{\mathit{b}}^{\mathit{n}}{\mathit{f}}^{\mathit{b}}+{\mathit{g}}^{\mathit{n}}+{\mathit{w}}_2$ ${\dot{\mathit{q}}}^{\mathit{n}}=\frac{1}{2}{\mathit{q}}^{\mathit{n}}{\omega }^{\mathit{b}},$ where

• Position in the navigation frame ${\mathit{t}}^{\mathit{n}}=[\mathit{x},\mathit{y},\mathit{z}]{}^{\mathrm{T}}$.

• Velocity in the navigation frame ${\mathit{v}}^{\mathit{n}}=[{\mathit{v}}_{\mathit{x}},{\mathit{v}}_{\mathit{y}},{\mathit{v}}_{\mathit{z}}]{}^{\mathrm{T}}$.

• Accelerometer measurement ${\mathit{f}}^{\mathit{b}}=[{\mathit{f}}_{\mathit{x}},{\mathit{f}}_{\mathit{y}},{\mathit{f}}_{\mathit{z}}]{}^{\mathrm{T}}$.

• ${\mathit{g}}^{\mathit{n}}$ is the gravitational acceleration.

• ${\mathit{R}}_{\mathit{b}}^{\mathit{n}}$ is a transforming matrix of a vector from the body to navigation frame.

• Attitude quaternion ${\mathit{q}}^{\mathit{n}}=[{\mathit{q}}_0,{\mathit{q}}_1,{\mathit{q}}_2,{\mathit{q}}_3]{}^{\mathrm{T}}$.

• Gyroscope measurement ${\omega }^{\mathit{b}}=[0,{\omega }_{\mathit{x}},{\omega }_{\mathit{y}},{\omega }_{\mathit{z}}]{}^{\mathrm{T}}$ (in a quaternion form).

• w ₁ and w ₂ are the white Gaussian process noises with strength matrices of Q ₁ and Q ₂ for the translation and velocity, respectively.

The runway pose solutions $\mathit{z}$ from the previous section are subsequently used as the observations to the extended Kalman filter, which are related to the inertial navigation state $\mathit{z}=\left[\begin{array}{c}\mathit{t}\\ \mathit{q}\end{array}\right]=\left[\begin{array}{c}{\mathit{R}}_{\mathit{b}}{\mathit{R}}_{\mathit{n}}^{\mathit{b}}{\mathit{t}}^{\mathit{n}}+{\mathit{v}}_1\\ \mathit{q}({\mathit{R}}_{\mathit{b}}{\mathit{R}}_{\mathit{n}}^{\mathit{b}})+{\mathit{v}}_2\end{array}\right],$ where v ₁ and v ₂ are the white Gaussian observation noises with strength matrices of R ₁ and R ₂ for the translation and quaternion, respectively. By using the non-linear process and observation models and the corresponding Jacobian matrices, the state vector and its covariance are propagated and updated (Fig. 3).

[IMAGE OMITTED. SEE PDF]

Experiment results

FlightGear interface

The FlightGear package is a continuously developing open-source flight simulator (). During runtime, rendered images were extracted from FlightGear and saved in an image format, while the instantaneous state information was stored as one large text. The images were also displayed in the open-source GUI display tool, EZWGL. This allowed the functionality of algorithms to be monitored during run time. A typical output is illustrated in Fig. 1. Two images are shown in the figure. On the left is the raw grey-scale image extracted from FlightGear and on the right is the same image with its corners highlighted. The four runway corners are clearly marked in addition to a multitude of white runway marking corners. Although only the corners marking the extremities of the runway were used in this application, the additional runway corners would serve to improve the state estimate. During the simulation, the aircraft was controlled by the user to make several lengthwise passes over the runway while ensuring the complete runway appeared in the image for as long as possible. Corner extraction was performed using the SUSAN [20] corner detecting algorithm, chosen for its high performance, adaptability and robustness to noise. For realism noise was inserted to body-derived inertial measurements. If the complete runway was contained in the image, then the set of camera-to-corner unit vectors was ascertained. The discretisation and rendering process of the simulator precluded the necessity for additional noise to be inserted on these measurements.

Results

Fig. 4 compares the estimated and actual 3D and 2D paths traversed by the aircraft showing the correctional effects when the runway was in full view of the camera. Fig. 5 (left) further compares the estimated aircraft xyz -position with the actual obtained during the experiment. The effect of adding Gaussian noise to the accelerometer and gyroscope data is most noticeable in the position plot. This is clearly evident in the aircraft's y -position, as the y -position appears to be affected the most because of the low magnitudes relative to the other two axes. Fig. 5 (right) illustrates the number of iterations required for convergence of the LM process through the usable range of images. Convergence was assumed when the residual was less than the value of ${10}^{-2}$. Thresholds lower than this are not guaranteed to produce better results. A weighting of 100 was placed on the $({\mathit{q}}^{\mathrm{T}}\mathit{q}-1)$ term in the cost function for the faster convergence. Fig. 6 shows the estimated errors in velocity and attitude (converted to Euler angles for display) with $1\sigma$ uncertainty. Fig. 7 provides the filter's innovation sequences of the position and attitude with $1\sigma$ uncertainty, showing the consistent operations of the filter.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Conclusions

This work presented the development of an algorithm that allows the use of visual runway observations as zero bias, absolute navigation references. Runways are suitable for this task because their geodetic coordinates are known precisely, as are their dimensions and detection markings, and the fact that they are designed to be highly visible. The emphasis in this application was placed on increasing the level of situational awareness for an aircraft during its landing approach, providing additional redundancy to navigating in the vicinity of a runway, thereby reducing the reliance on global positioning system. The algorithm produces a UDQ-based cost function based on the geometric relationships between the observations and the runway, and an iterative minimisation of the cost function using the LM technique. The algorithm was then integrated into the inertial navigation system and demonstrated using a flight simulator using a camera-in-the-loop setup. Future work is extending the method to accommodate the state-of-the-art non-linear observer-based SLAM technique and the UDQ to enhance stability and efficiency.

Acknowledgments

The work presented here is an extract from the primary author's PhD thesis, Vision-Aided Aircraft Navigation and Trajectory Optimisation, completed at the Australian National University in 2007.

References

1

Kim J. Byun H. Guivant: ‘Compressed Pseudo‐SLAM: Pseudorange Integrated Generalised Compressed SLAM’, Australasian Conference of Robotics and Automation, Brisbane, December 2020

2

Usenko V. Demmel N. Schubert D. et al.: ‘Visual‐inertial mapping with non‐linear factor recovery’, IEEE Robot. Autom. Lett., 2020, 5, (2), pp. 422–429

3

Liang Q. Liu M.: ‘A Tightly Coupled VLC‐Inertial Localization System by EKF’, IEEE Robot. Autom. Lett., 2020, 5, (2), pp. 3129–3136

4

Forster C. Carlone L. Dellaert F. et al.: ‘On‐manifold preintegration for real‐time visual–inertial odometry’, IEEE Trans. Robot., 2017, 33, (1), pp. 1–21

5

Miiller M.G. Steidle F. Schuster M.J. et al.: ‘Robust visual‐inertial state estimation with multiple odometries and efficient mapping on an MAV with ultrawide FOV stereo vision’. 2018 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS), Madrid, Spain, October 2018, pp. 3701–3708

6

Ashford R.: ‘Flight safety digest: a study of fatal approach‐andlanding accidents worldwide 1980–1996’, Flight Safety Digest, 1998, 17, (2/3)

7

Serra P. Cunha R. Hamel T. et al.: ‘Landing of a quadrotor on a moving target using dynamic image based visual servo control’, IEEE Trans. Robot., 2016, 32, (6), pp. 1524–1535

8

Xu L. Luo H.: ‘Towards autonomous tracking and landing on moving target’. 2016 IEEE Int. Conf. on Real‐time Computing and Robotics (RCAR), Angkor Wat, Cambodia, June 2016, pp. 620–628

9

Tang D. Jiao Y. Chen J.: ‘On automatic landing system for carrier plane based on integration of INS, GPS and vision’. 2016 IEEE Chinese Guidance, Navigation and Control Conf. (CGNCC), Nanjing, China, August 2016, pp. 2260–2264

10

Nguyen K.D. Ha C.: ‘Vision‐based hardware‐in‐the‐loop simulation for unmanned aerial vehicles’, in ‘Intelligent computing theories and application’ (Springer International Publishing, Wuhan, China, 2018), pp. 72–83

11

Coutard L. Chaumette F. Pflimlin J.‐M.: ‘Automatic landing on aircraft carrier by visual servoing’. 2011 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, San Francisco, USA, 2011, pp. 2843–2848

12

Kus M.: ‘Autonomous carrier landing of a fixed‐wing UAV with airborne deck motion estimation’. Ph.D. dissertation, 2019

13

Gibert V. Plestan F. Burlion L. et al.: ‘Visual estimation of deviations for the civil aircraft landing’, Control Eng. Pract., 2018, 75, pp. 17–25

14

Moore A.J. Schubert M. Dolph C. et al.: ‘Machine vision identification of airport runways with visible and infrared videos’, J. Aerosp. Inf. Syst., 2016, 13, pp. 266–277

15

Wu Y. Hu X. Hu D. et al.: ‘Strapdown inertial navigation system algorithms based on dual quaternions’, IEEE Trans. Aerosp. Electron. Syst., 2005, 41, (1), pp. 110–132

16

Wang X. Yu C. Lin Z.: ‘A dual quaternion solution to attitude and position control for rigid‐body coordination’, IEEE Trans. Robot., 2012, 28, (5), pp. 1162–1170

17

Cheng J. Kim J. Shao J. et al.: ‘Robust linear pose graph‐based SLAM’, Robot. Auton. Syst., 2015, 72, pp. 71–82

18

Han P. Cheng Z. Chang L.: ‘Automatic runway detection based on unsupervised classification in polsar image’. 2016 Integrated Communications Navigation and Surveillance (ICNS), Herndon, USA, April 2016, pp. 6E3‐1–6E3‐8

19

Abu‐Jbara K. Alheadary W. Sundaramorthi G. et al.: ‘A robust vision‐based runway detection and tracking algorithm for automatic UAV landing’. 2015 Int. Conf. on unmanned aircraft systems (ICUAS), Denver, USA, 2015, pp. 1148–1157

20

Smith S.M. Brady J.M.: ‘SUSAN ‐ a new approach to low level image processing’. Int. Journal of Computer Vision, 1997, 23, (1), pp. 45–78