**Alexander Miller, Boris Miller \*, Alexey Popov and Karen Stepanyan**

Institute for Information Transmission Problems RAS, Bolshoy Karetny per. 19, Build.1, Moscow 127051, Russia; amiller@iitp.ru (A.M.); ap@iitp.ru (A.P.); KVStepanyan@iitp.ru (K.S.)

**\*** Correspondence: bmiller@iitp.ru; Tel.: +7-495-6504781; Fax: +7-495-6500579

Received: 3 January 2019; Accepted: 10 March 2019; Published: 18 March 2019

**Abstract:** An automatic landing of an unmanned aerial vehicle (UAV) is a non-trivial task requiring a solution of a variety of technical and computational problems. The most important is the precise determination of altitude, especially at the final stage of approaching to the earth. With current altimeters, the magnitude of measurement errors at the final phase of the descent may be unacceptably high for constructing an algorithm for controlling the landing manoeuvre. Therefore, it is desirable to have an additional sensor, which makes possible to estimate the height above the surface of the runway. It is possible to estimate all linear and angular UAV velocities simultaneously with the help of so-called optical flow (OF), determined by the sequence of images recorded by an onboard camera, however in pixel scale. To transform them into the real metrical values it is necessary to know the current flight altitude and the camera angular position values. The critical feature of the OF is its susceptibility to the camera resolution and the shift rate of the observed scene. During the descent phase of flight, these parameters change at least one hundred times together with the altitude. Therefore, for reliable application of the OF one needs to coordinate the shooting parameters with the current altitude. However, in case of the altimeter fault presence, the altitude is also still to be estimated with the aid of the OF, so one needs to have another tool for the camera control. One of the possible and straightforward ways is the camera resolution change by pixels averaging in computer part which performed in coordination with theoretically estimated and measured OF velocity. The article presents results of such algorithms testing from real video sequences obtained in flights with different approaches to the runway with simultaneous recording of telemetry and video data.

**Keywords:** UAV; landing; optical flow; video navigation; Kalman filter

### **1. Introduction**

Today the development of optoelectronic devices and data transmission systems allows to use them in remote control systems of unmanned aerial vehicles (UAV). In case of remote control which is carried out by an operator the characteristics of the optical devices and the transmitted image must correspond to the capabilities of human vision. However, in case of the UAV autonomous flight, the optical systems and the onboard computer must work together solving the problems of identifying the observed objects with their coordinates. It means that requirements to characteristics of optoelectronic devices are different from remote flight control case. The UAV control system which includes an optoelectronic system (OES) and an onboard computer determines the movement of the onboard video camera and identifies the objects observed in the field of view of the camera [1]. There are several approaches to the usage of OES [2]. The first one is to detect and to track the movement of specific local areas (reference points) on the image by analogy with human vision [3,4]. With this approach, it is easy to transform the recorded images into metric values for the control system. Some examples of this approach to UAV navigation are as follows:


One possible reference to this approach is *sparse OF*. An example of this approach to the UAV landing in case of the smoke occlusions basing on the natural terrain landmarks is given in [11]. The second approach is based on non-metric analysis and is built by analogy with a vision of insects or birds [12]. Nowadays an approach that uses information containing in the vector field of the image motion (in other words in the optical flow (OF)) is developing. OF is well known to engineers developing cameras for shooting from moving carriers where unreduced OF leads to the resolution degradation and needs reduction either optomechanically [13,14] or electronically [15]. At the same time in the UAV application area, the OF contains the information about the carrier's velocities and thereby may serve as an additional sensor for the UAV navigation.

OF is the projection of the camera's motion onto the focal plane. OF generates the nonuniform field of image shifts velocities which is known as *dense optical flow*. In case of tracking the displacement of some reference points on the underlying surface, it is common to use the *sparse optical flow* term. The methodology of the OF computation had been developed long ago, mainly for estimation of quality for various optomechanical image shift compensation systems performance [16]. Nowadays there are various examples of the OF usage in UAV applications such as:


Our research team is working on video navigation as an additional navigation tool for UAV of aeroplane type. We consider OF as an additional tool for estimation of the UAV velocity in the absence of specific regions on the earth surface which can serve as beacons for estimation of the UAV position. The specific feature of the OF is that it gives information about velocities only but not on position. Therefore, the bias of the position estimates which inevitably exists at the beginning of the flight path can only increase without additional intermediate corrections. Correct filtering algorithm can reduce this bias but cannot eradicate it. It means that it is essential to evaluate the specific noises related to the OF estimation, so in our experiments, we perform flights where qualified pilot performs the series of approaches to the runway, during which one can carefully record the video sequences and corresponding telemetry data. These data serve as a source of information about specific OF noises. Filtering method fuses the data at UAV altitude estimation during landing.

In the existing literature, the OF usage at landing generally relates to the copter type UAV [11]. Here it is possible to coordinate velocity of descent with the observable video sequence. There are a series of works demonstrating successful approaches to the copter control based on the divergence observation. Here the divergence of the OF field serves as a measure of the approach velocity to the earth surface. Almost all articles were presenting the successful application of the OF relate to micro air vehicles (MAV), where OF used with supplemental range meter. For example, in [26] authors consider a Vertical-Take-OFF-and-Landing UAV (VTOL UAV) with an Inertial Measurement Unit

(IMU) equipped with suitable filtering algorithms providing reliable estimates of the UAV pose and rotational velocities. Moreover, this UAV equipped with an optoelectronic camera serving as an additional passive sensor whose output is enriching in-flight information. These authors also assume that the UAV control system has an image processing capability to identify a target and to compute OF over the full image basing on well-known Horn and Schunck approach [27].

In [28] a nonlinear controller is presented for a VTOL UAV that exploits a measurement of the average OF to enable hovering and landing onto a moving platform such as the deck of a seagoing vessel.

In a detailed review of copter type UAV's [11] video control presents a variety of approaches to the usage of natural landmarks and the OF approaches. The principal differences between UAVs of a copter type and even MAV and standard sized UAVs lay in their sizes and velocities. When small-sized UAV is approaching the obstacles with low velocity, the latter can be controlled by the OF signal, since the divergence of the OF field serves as a measure of approaching speed [29–32]. On the contrary, the standard sized UAVs are landing onto the runway and approaching to the earth by standard glissade during which the altitude changes from hundreds of meters to zero and velocity reduces from dozens of meters per second to meters per second and finally to zero. Simultaneously the rate of the image motion changes on the same orders which lead to the resolution degradation, so the detection of the natural landmarks like in [11] becomes impossible. An application of the OF techniques which determine the image velocity via analysis of the evolution of local image areas needs the coordination of the resolution level of OES with the current speed. The measure of such coordination is obtainable by comparison of measured and theoretically calculated OF velocities.

Even if there are communications related to successful usage of the OF with the aid of two vertically displaced cameras [33,34], more careful analysis shows that with real images and flights the bias between the estimated and real values of the OF achieves unacceptable values [35].

That is why the approach developed for MAV do not apply to the aeroplane landing where the descent velocity cannot be arbitrary controlled. Moreover, the OES resolution when approaching to the earth changes about a hundred times from the start of a glissade. By this reason, the standard OF evaluation does not work accurately if the OES resolution does not reduce in coordination with the UAV altitude. A changing of the camera focal length is somewhat tricky since it needs an additional controlled optomechanical system. Meanwhile, the resolution can be changed by averaging pixels in coordination with estimated OF velocity obtained from the video sequence and its comparison with the value calculated theoretically from exact OF formulas based on current altitude estimation. The principal aim of this article is to demonstrate this approach via real video sequences obtained in real flights. There are two ways to use video navigation in autonomous UAV flights. The first one is a detection of the terrain objects with known coordinates, a definition of their aspect angles and finally, a determining the current UAV coordinates. We demonstrated this approach in [5]. However, if such objects are absent in the field of view, one can determine the current velocity and coordinates by filtering, and the OF is a good way for that. So the combination of direct measurements and OF provides the continuous tracking of the UAV trajectory. This way is convenient in so-called cruise part of the flight, where the altitude and velocity remain almost constant. The landing is the different issue, since the altitude changes on few orders as well as the scale of the picture, therefore, small details become apparent and affect the accuracy of the determining of the OF [36]. That is why at landing phase one needs to adapt the camera resolution and the frame rate following the current altitude. Creation of such cameras is a separate problem, though one can change the resolution by averaging the image signal. In this article, we show how to create such averaging algorithm, which together with filtering based estimation guarantees the reliable evaluation of the altitude at the descent from the altitude of 300 m to 5 m.

The structure of the article is as follows. In Section 2, we present the theory of the OF computation. In Section 3, we give the algorithm of the UAV motion estimation with the aid of the OF and Kalman filtering. Section 4 gives the results of the OF obtained during the real UAV flight at the landing

phase and presents the approach to the resolution adaptation from current altitude estimation. In Section 5 we give the algorithms of the averaging scale switching based on the comparison of exact OF values (theoretically calculated) and their estimates with the aid of Lucas-Kanade (L–K) algorithm [37]. It shows that the difference between estimated and calculated OF values may serve as a sensor for the control of averaging. The results of numerical tests showing the reliable tracking of the altitude up to the 5m provided. Section 6 is conclusions where we discuss the results.

#### **2. OF Computation: Theory**

Since the end of the 70s, many mentions of the OF formulas appeared in the literature, see for example [38]. However, these formulas usually relate to the case of fixed camera position and cannot take into account possible inclination of the line of sight which occurs either during the plane turn (azimuth and roll angles) or during the landing-descent (pitch angle). Meanwhile, the case of a camera with a stabilised line of sight is much more complicated. As an example, the copter needs some inclination angle in order to create propulsion in the flight direction. Hence the pitch angle correction is necessary as it presented in the example in [39]. In recent years navigation by computation of the camera path and the distance to obstacles with the aid of field of image motion velocities (i.e., OF) became highly demanded particularly in the area of relatively small and even micro UAV. Video sequences captured by onboard camera give the possibility of the onboard OF calculation with the aid of relatively simple algorithms like Lucas-Kanade [37] and even more sophisticated ones using higher order terrain illuminance approximations [40–42].

Theory of the OF computation in general flight conditions is in [43–45]. On that basis, the specialised software IMODEL (version 8.2.o, proprietary) developed and various scenarios of the OF usage in the UAV navigation presented. The IMODEL permits to analyse both the *direct* and *inverse* problems of the OF computation. The *direct* problem is the calculation of the image motion at any point of the field of view for the general orientation of the line of sight. The *inverse* one is the estimation of the OF field for given modelled moving landscape registered by the camera virtually. Moreover, in the implementation stage, the calculus precision problem arises. On one side the higher level filtering algorithm described in [45] and the exact parameters for it affect the resulting precision of the estimated parameters of UAV position. On the other side, the precision achieved is the consequence of the chosen low-level method to determine OF.

An example of the software application is the evaluation of the altitude estimation algorithm with the aid of two vertically displaced cameras [33,35] which shows the presence of bias in the altitude estimation.

A general approach to the OF computation is based on the pinhole camera geometrical model assuming the flight over a flat surface [39,43,44] (see Figures 1 and 2).

**Figure 1.** The picture shows the projection of the image plane onto the earth surface. The orientation of the line of sight is changing by rotation at point G.

We use the camera coordinate system (*ξ*, *η*, −*F*), where (*ξ*, *η*) are the pixel coordinates, *ξ* corresponds to the direction of flight, *η* is perpendicular to *ξ*, and *F* is the lens focal length, and the third coordinate directs to the principal camera point. As it follows from the general theory [16,45] the OF velocities (*V<sup>ξ</sup>* , *Vη*) equal to

$$
\begin{pmatrix} V\_{\xi} \\\\ V\_{\eta} \end{pmatrix} = \begin{pmatrix} \frac{\partial \tilde{\xi}}{\partial t} \\\\ \frac{\partial \eta}{\partial t} \end{pmatrix} \bigg|\_{\mathbf{x} = \mathbf{x}(\xi, \eta, \Lambda(t)) : \mathcal{Y} = \mathcal{Y}(\xi, \eta, \Lambda(t))} \tag{1}
$$

where dependence on Λ(*t*) includes the current attitude of the vehicle and orientation of the line of sight. For in-flight computations, one needs to use more extended formula including *Wx*, *Wy*, *WH*, *ωp*, *ω<sup>r</sup>* , *ω<sup>y</sup>* which is the following matrix equation:

$$
\begin{pmatrix} V\_{\tilde{\xi}} \\\\ V\_{\eta} \end{pmatrix} = D\_1(\tilde{\xi}, \eta, \Lambda(t)) \begin{pmatrix} W\_{\mathcal{X}} \\\\ W\_{\mathcal{Y}} \\\\ W\_H \end{pmatrix} + D\_2(\tilde{\xi}, \eta, \Lambda(t)) \begin{pmatrix} \omega\_p \\\\ \omega\_r \\\\ \omega\_y \end{pmatrix} \tag{2}
$$

where (*Wx*, *Wy*) *<sup>T</sup>* are the velocities components of the UAV in the plane parallel to the earth surface, (*ωp*, *ω<sup>r</sup>* , *ωy*) *<sup>T</sup>* are the rotation velocities relative to the line of sight, and *W<sup>H</sup>* is the vertical velocity component of the UAV. Exact formulas may found in [39]. By using the estimations of the left-hand side (LHS) of (2) either via Lucas-Kanade (L–K) [37] algorithm or any other, one gets the additional sensor for the UAV velocities estimation. The principal difficulties are the implicit dependence on Λ(*t*) including the line of sight orientation, given by angles (*p*,*r*, *y*) conditionally called the *pitch, roll, yaw*, and the current altitude *H*.

Therefore, the observation algorithm requires estimation of all mentioned above parameters by Kalman filtering via observation of their derivatives from (2).

The system (2) contains observable variables in the LHS and a set of six variables to be estimated in the right-hand side (RHS). For a sufficiently large field of view, one can evaluate velocities in a vast amount of pixels, while the values are the same for any point. It gives the possibility to solve the issue for example via least square method since the observation model is linear in velocities, and it is possible to estimate unobservable variables via Kalman filtering [45].

**Remark 1.** *In our approach, we use OF just as an additional sensor to estimate the altitude only. The reason is that in these experiments we use the data obtained with the aid of a series of landing approaches made by a qualified pilot of the light aeroplane. Other parameters such as velocity and orientation angles can be obtained from inertial navigation system (INS) with sufficient accuracy, so we use their nominal values with small corrections made with the aid of Kalman filter (see below* (5)*,* (6)*). The current value of the altitude is necessary for estimation via OF. In the real flight of the UAV all possible navigation sensors must be used in the fusion with the OF. Here we just test the ability of the OF solely, therefore the real approaches to the runway have been performed by a pilot of a light aeroplane with recording of the video images synchronized with telemetry data from INS and satellite navigation system (SNS).*

#### **3. Estimation of the UAV Motion by the OF and the Kalman Filtering**

The model of UAV motion described below. It includes the UAV dynamic model and generic measurements model which based on the OF estimation of the UAV attitude velocities.

#### *3.1. The UAV Linear Velocity Estimation*

The UAV velocity vector **V** = *col*(*Vx*, *Vy*, *Vz*) by coordinates *x*, *y*, *z* :

$$\mathbf{V}(t\_{k+1}) = \mathbf{V}(t\_k) + \mathbf{a}(t\_k)\Delta t + \mathbf{W}(t\_k) \tag{3}$$

where *t<sup>k</sup>* is the current time, *t<sup>k</sup>* = *t*<sup>0</sup> + *k*∆*t*,

> **a**(*t<sup>k</sup>* ) = *col*(*ax*, *ay*, *az*) — the vector of accelerations coming from INS (inertial navigation system), **W**(*t<sup>k</sup>* ) — is the vector of the current perturbations in UAV motion.

We assume that the components of the perturbation vector are white noises with variances (*σ* 2 *x* , *σ* 2 *y* , *σ* 2 *z* ).

Velocity measurements using OF have the following general form:

$$\mathbf{m}\_V(t\_k) = \mathbf{V}(t\_k) + \mathbf{W}\_V(t\_k),\tag{4}$$

),

where **W***V*(*t<sup>k</sup>* ) — are uncorrelated white noises with variances (*σ* 2 *Vx* , *σ* 2 *Vy* , *σ* 2 *Vz* ). Consider relations (3) and (4) for the velocity along *x* axis:

$$\begin{aligned} V\_{\mathbf{x}}(t\_{k+1}) &= V\_{\mathbf{x}}(t\_k) + a\_{\mathbf{x}}(t\_k)\Delta t + \mathcal{W}\_{\mathbf{x}}(t\_k), \\\\ \mathcal{m}\_{V\_{\mathbf{x}}}(t\_k) &= V\_{\mathbf{x}}(t\_k) + \mathcal{W}\_{V\_{\mathbf{x}}}(t\_k). \end{aligned}$$

Velocity along *x* axis estimation on the *k* + 1 step:

$$\begin{aligned} \hat{\mathcal{V}}\_{\mathbf{x}}(t\_{k+1}) &= K\_{\mathbf{x}}(t\_{k+1}) m\_{V\_{\mathbf{x}}}(t\_{k+1}) + (1 - K\_{\mathbf{x}}(t\_{k+1})) \tilde{\mathcal{V}}\_{\mathbf{x}}(t\_{k+1}), \\\\ \hat{\mathcal{V}}\_{\mathbf{x}}(t\_{k+1}) &= \hat{\mathcal{V}}\_{\mathbf{x}}(t\_{k}) + a\_{\mathbf{x}}(t\_{k}) \Delta t. \end{aligned}$$

Kalman filter gives the estimation of *V*ˆ *x*:

$$\begin{split} \dot{\mathcal{V}}\_{\mathbf{x}}(t\_{k+1}) &= \mathcal{K}\_{\mathbf{x}}(t\_{k+1})m\_{V\_{\mathbf{x}}}(t\_{k+1}) + (1 - \mathcal{K}\_{\mathbf{x}}(t\_{k+1}))(\dot{\mathcal{V}}\_{\mathbf{x}}(t\_{k}) + a\_{\mathbf{x}}(t\_{k})\Delta t), \\ \mathcal{K}\_{\mathbf{x}}(t\_{k+1}) &= \frac{\dot{\mathcal{P}}^{V\_{\mathbf{x}}V\_{\mathbf{x}}}(t\_{k}) + \sigma\_{\mathbf{x}}^{2}}{\dot{\mathcal{P}}^{V\_{\mathbf{x}}V\_{\mathbf{x}}}(t\_{k}) + \sigma\_{\mathbf{x}}^{2} + \sigma\_{V\_{\mathbf{x}}}^{2}}, \\ \mathcal{P}^{V\_{\mathbf{x}}V\_{\mathbf{x}}}(t\_{k+1}) &= \frac{\sigma\_{V\_{\mathbf{x}}}^{2}\left(\dot{\mathcal{P}}^{V\_{\mathbf{x}}V\_{\mathbf{x}}}(t\_{k}) + \sigma\_{\mathbf{x}}^{2}\right)}{\dot{\mathcal{P}}^{V\_{\mathbf{x}}V\_{\mathbf{x}}}(t\_{k}) + \sigma\_{\mathbf{x}}^{2} + \sigma\_{V\_{\mathbf{x}}}^{2}}. \end{split} \tag{5}$$

The formulae for *V*ˆ *<sup>y</sup>* and *V*ˆ *<sup>z</sup>* are analogous.

**Remark 2.** *Kalman filter coefficients have been derived from empirically registered standard deviations of noise processes. They assumed to be constant during the whole duration of landing. Only the scale rate changes, but after application of the least squares method for estimation of velocities in the RHS of* (1) *for the large number of observations the remaining noise in observation of the LHS of* (1) *could be assumed constant for various averaging scales.*

#### *3.2. The UAV Angles and Angular Velocities Estimation*

UAV angular position estimation is given by three angles *θ*(*t<sup>k</sup>* ), *ϕ*(*t<sup>k</sup>* ), *γ*(*t<sup>k</sup>* ) that are (pitch, roll and yaw, respectively), angular velocities *ωp*(*t<sup>k</sup>* ), *ωr*(*t<sup>k</sup>* ), *ωy*(*t<sup>k</sup>* ) and angular accelerations *ap*(*t<sup>k</sup>* ), *ar*(*t<sup>k</sup>* ), *ay*(*t<sup>k</sup>* ).

Pitch angle and pitch angular velocity dynamics described by the following relations:

$$
\theta(t\_{k+1}) = \theta(t\_k) + \omega\_p(t\_k)\Delta t + a\_p(t\_k)\frac{\Delta t^2}{2}.
$$

$$
\omega\_p(t\_{k+1}) = \omega\_p(t\_k) + a\_p(t\_k)\Delta t + \mathcal{W}\_p(t\_k).
$$

where *Wp*(*t<sup>k</sup>* ) — is the white noise with variance *σ* 2 *p* .

The pitch angular velocity measurement using the OF has the following form:

$$
\boldsymbol{m}\_p(t\_k) = \boldsymbol{\omega}\_p(t\_k) + \mathcal{W}\_{\boldsymbol{\omega}\_p}(t\_k)\_\prime
$$

where *Wω<sup>p</sup>* (*tk* ) — is the noise in the angular velocity measurements using OF, which is the white noise with variance *σ* 2 *ωp* .

Similarly to the linear velocity estimation we get the pitch angle *θ*(*t<sup>k</sup>* ) and pitch angular velocity *ωp*(*t<sup>k</sup>* ) estimations:

$$\begin{split} \dot{\theta}(t\_{k+1}) &= \dot{\theta}(t\_k) + \hat{\omega}\_p(t\_k)\Delta t + a\_p(t\_k)\frac{\Delta t^2}{2} \\\\ \hat{\omega}\_p(t\_{k+1}) &= K\_p(t\_{k+1})m\_p(t\_{k+1}) + (1 - K\_p(t\_{k+1}))(\hat{\omega}\_p(t\_k) + a\_p(t\_k)\Delta t), \\\\ K\_p(t\_{k+1}) &= \frac{\hat{P}^{\omega\_p\omega\_p}(t\_k) + \sigma\_p^2}{\hat{P}^{\omega\_p\omega\_p}(t\_k) + \sigma\_p^2 + \sigma\_{\omega\_p}^2}, \\\\ \hat{P}^{\omega\_p\omega\_p}(t\_{k+1}) &= \frac{\sigma\_{\omega\_p}^2(\hat{P}^{\omega\_p\omega\_p}(t\_k) + \sigma\_p^2)}{\hat{P}^{\omega\_p\omega\_p}(t\_k) + \sigma\_p^2 + \sigma\_{\omega\_p}^2}. \end{split} \tag{6}$$

The formulae for *ϕ*ˆ, *γ*ˆ and *ω*ˆ*<sup>r</sup>* , *ω*ˆ*<sup>y</sup>* are analogous.

#### *3.3. Joint Estimation of the UAV Attitude*

Jointly (5) and (6) give the estimation **Vˆ** of the attitude parameters vector, namely:

$$\mathbf{V} = \operatorname{col}(V\_{\mathbf{x}\prime}V\_{y\prime}V\_{z\prime}\boldsymbol{\omega}\_{p\prime}\boldsymbol{\omega}\_{r\prime}\boldsymbol{\omega}\_{y})\dots$$

This vector measured via OF for each pixel in each frame, which gives according to (2) the overdetermined system of linear equation for **V** entries.

Since all noises are uncorrelated the covariance matrix *P* for **V** is diagonal with entries (5), (6).

Integration of velocities **V** gives the estimation of the UAV position which used for estimation of matrices *D*1, *D*<sup>2</sup> in (2).

#### *3.4. Discussion of the Algorithm*

Many uncertainties remain in the algorithm described above. They may be evaluated only in the real flight with real data. The principal uncertainty is the noise level in formulas (2) and its dependence on the altitude of flight. Experiments show that when approaching to the earth, the tiny details become very important and corrupt the OF estimation. In order to suppress these details, one needs to change resolution but to keep velocity estimation to be correct. In the real flight by averaging the pixel samples one can change resolution, but it is necessary to do it in coordination with the current altitude. Our experiments show that it is possible to change the resolution by comparison of the L–K OF estimation and theoretically calculated OF velocity through the current estimation of the altitude. That is a heuristic approach, and only experiments with real data may show either is it appropriate or not. The difference between two ways of the OF calculation is that L–K does not use the data from INS and even the estimation of the altitude obtained on previous steps with the aid of Kalman filtering. On the other hand, the OF calculation from theoretical model uses the current UAV attitude estimation. In the next section, we demonstrate an approach to the development of such an algorithm to control the scaling.

### **4. OF Estimation in the Real Flight**
