Prediction of microunmanned aerial vehicle flight behavior from two-dimensional intensity images

Abstract. The increasing number of microunmanned aerial vehicles (MUAVs) is a rising risk for personal privacy and security of sensitive areas. Owing to the highly agile maneuverability and small cross section of the MUAV, effective countermeasures (CMs) are hard to deploy, especially when a certain temporal delay occurs between the localization and the CM effect. Here, a reliable prediction of the MUAV flight behavior can increase the effectiveness of CMs. We propose a pose estimation approach to derive the three-dimensional (3-D) flight path from a stream of two-dimensional intensity images. The pose estimation in a single image results in an estimation of the current position and orientation of the quadcopter in 3-D space. Combined with flight behavior model, this information is used to reconstruct the flight path and to predict the flight behavior of the MUAV. In our laboratory experiments, we obtained a standard deviation between 1 and 24 cm in a five-frame prediction of the 3-D position, depending in the actual flight behavior.


Introduction
In the recent past, the number of incidents with intentional and unintentional misuse of small or microunmanned aerial vehicles (MUAVs) has increased due to the increasing number of commercially available small quadcopters.In their annual reports, 1 the Federal Aviation Administration has pointed out an increasing number of sightings and close encounters of MUAVs reported in civilian air traffic.Further, the first serious incidents where a MUAV crashed into an air plane have occured. 2Beside this threat to civilian air traffic, misuse of MUAVs can also harm personal and public privacy and create security issues.
4][5][6][7] Thus, in the past decade, several groups have been working on the problem of efficiently detecting, tracking, and classifying MUAVs in different scenarios.5][16] Classification can also be used to distinguish MUAVs from, e.g., birds.
Different sensors (such as RADAR, acoustics, and optics) will be embedded in a counter-MUAV system where sensor information will be processed to reliably detect/track, classify, and localize MUAV.Based on this information, decisions are made and countermeasures (CMs) are engaged.Owing to the internal processing time of information generation and decision-making as well as external time constrains (e.g., propagation of a projectile), the effect on the MUAV will occur after a certain inherent time delay Δt.Owing to the agile operation and small cross section, the probability of successful CMs is relatively low. 17Here, the prediction of the MUAV position at the moment of effect could significantly increase the success rate.
First, works that predict the MUAV position in the image plane of the succeeding frame (frame-to-frame prediction) of a video stream 18,19 using a linear motion model have been presented.In our approach, we predict the position of the MUAV in three-dimensional (3-D) space and on a longer time scale t P (with t P ≈ Δt) or over n-frames.For a reliable result, we determine the 3-D position of the MUAV from pose estimation in a two-dimensional (2-D) intensity images and apply a navigation model to simulate the flight behavior.

Prediction of Microunmanned Aerial Vehicle
Position in the Image Plane Based on linear quadratic estimation, Kálmán 20 developed a filter to either reduce measurement noise effects and to estimate values (or determine the most probable value).In the computer vision community, for instance, the Kalman filter approach is widely used for tracking and navigating the MUAV. 18,19The Kalman filter fð•Þ is used to predict the MUAV state vector X i , which contains a set of parameters, from a current measurement y i and the previous state vector X i−1 , as defined in Eq. (1).Here, ω describes noise characteristics E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 1 6 4 X i ¼ fðX i−1 ; y i ; ωÞ: (1) Depending on the application (e.g., navigation or tracking), the state vector summarizes several physical parameters [such as position, velocity, angles (roll, yaw, pitch), and angular velocities] to describe the current flight state.In our approach, the state vector is defined as X i ¼ ½x i ; ẋi ; ẍi T containing only the MUAV position x i , velocity _ x i , and acceleration ẍi that are the first and second temporal derivatives of the position in the i'th frame.Higher order of temporal derivatives [the jerk (third) and the snap (fourth)] could be of interest in future algorithms and are neglected in the current code.
We developed a tracking and prediction filter based on Murphy's "Kalman filter toolbox for Matlab" 21 and the tutorial of Welsh and Bishop. 22Our filter was designed to predict a MUAV position in a n-frame future from a continuous data stream and has two stages: a Kalman filter to estimate the current state vector X i and an n-frame prediction filter to estimate the future state Xiþn .
The Kalman filter, as described in Sec. 5, returns a state vector X i at time t i , containing position, velocity, and acceleration.Then, the MUAV state vector Xiþn at time t iþn can be predicted using the motion equation approach, as defined in Eq. ( 2).The predicted position xiþn can be extracted by Eq. ( 3).We define the motion prediction equations as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 5 5 1 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 4 9 3 The complete n-frame prediction procedure is illustrated in Algorithm 1 as a pseudocode.The computation process returns a predicted position ẑiþn and needs a set of input data consisting of the current image frame F i , the previous MUAV state vector X i−1 , and the number of frames n as the prediction range.First, the MUAV position x i in the current image has to be determined by, for instance, a tracking algorithm.This step will be discussed in the next paragraph.From these position measurements, the current state vector is calculated using the Kalman filter and the MUAV state Xiþn at frame i þ n is predicted.Finally, the estimated future position xiþn is extracted [see Eqs. ( 2) and ( 3)].
In Fig. 1, some frames of a video stream are presented.This video was recorded inside our laboratory, showing a small quadcopter MUAV (type DJI Phantom 3) flying around in the camera field of view.The MUAV is tracked with a state-of-the-art tracker based on the consensusbased matching and tracking of key points (CMT). 23,24n this algorithm, the image is transformed to an alternative base using feature detection from accelerated segment test (FAST) 25,26 and binary robust invariant scalable key points of features 27 as feature detector and descriptor subroutines.
The tracker employs clustering of correspondences as the central idea to distinguish inlier and outlier key points.This approach improves the tracking results, even when correspondences are located on deformed parts of the object.As the flying MUAV shows very rapid aspect changes due to the angle between the vision system and the target, this tracker is well adapted to this particular situation.The tracking algorithm returns a vector containing the MUAV position (x i and y i are the column and row positions of the MUAV, respectively) and the size of the bounding box (width and height: w i , h i ).In Fig. 1, some frames of the video stream illustrate the CMT tracking results.
An exemplary result of our tracking and prediction algorithm (Algorithm 1) is presented in Fig. 2. In Fig. 2(a), a fiveframe prediction result from our approach using first and second derivatives (□) is compared to a simple prediction approach using the first derivative only (Δ).Both prediction paths are compared to the ground-truth CMT tracking data (-).As long as the MUAV flight operation is slow and smooth, only small differences between ground truth and predicted paths are observed.But at turning points and during high-agile flight operation, the predictions are less accurate, especially for the simple prediction approach.
A deeper discussion of the prediction results is presented in Fig. 2(b).In the top diagram, the amount of single flight behavior components are plotted for each frame as the amplitude of the first and second derivatives [nj_ xj (⋯) and 1 2 n 2 jẍj (-)].It is obvious that rapid changes in the first derivative lead to peaks in the amplitude of the second derivative.
Furthermore, in the bottom diagram of Fig. 2(b), the prediction error is illustrated as the mismatch distance between ground truth and predicted position given by the L2-norm distance: d L2 ¼ kx i − xi k 2 .For the simple prediction approach (⋯), high errors can be observed in frames that are related to high agile or dynamic MUAV flight operations (see frames 110 to 125, 165 to 195, 205 to 220, etc.).In this specific sequence, we obtain a mean prediction error of 10.1 pixels.This mismatch distance can be reduced significantly by additionally employing the second derivative (-) in our prediction approach.With this algorithm it is possible to reduce all previously occurring peaks significantly and a mean prediction mismatch distance of 4.7 pixels is reached.9][30][31][32][33] These algorithms analyze the dynamic behavior of the spatial and the angular motion to control the power consumption of the four rotors individually.Typically, these models combine a global or inertial and a MUAV or body coordinate system.In our approach, we rely on the global coordinate system of the optical sensing device.In the coordinate system the x-and y-directions are the width and the height corresponding to the previous pixel row and column directions, respectively.Then, the z-direction is the range.
Further, we developed a simplified flight behavior model to describe the MUAV motion from the thrust T, the gravitational force G, and the aerodynamic drag force D, which can be derived from some fundamental parameters such as position and orientation in 3-D space, velocity, acceleration, and few boundary conditions (aerodynamic constants and parameters).In principle, we assume that the MUAV is a floating object (the vertical component of the thrust compensates for earth gravitation).Details of this fight behavior model can be found in Secs.6 and 7.
Position and orientation of the MUAV in 3-D space can be derived from pose estimation.For this task, the quadcopter used in our experiment can be modeled by a skeleton model consisting of four small vertical segments equally distributed around a bigger vertical segment, as shown in Fig. 3.The four small segments represent the rotation axes of each propeller and the bigger one represents the body of the quadcopter.
The use of a skeleton model has the benefit of not being too specific to the investigated MUAV.Later adaption to another MUAV can easily be done by changing the model.The only values that have to be known are the number of propellers and their distance to the center.
For the investigated vehicle, the DJI Phantom 3, 34 there are four propellers located 175 mm from the center.The number of points per segment can also vary, but it does not have a great influence; we usually use 20 points per segment.The skeleton points are designated by S ¼ fs j g when placed at the origin.The pose of the skeleton corresponding to the image i is represented by a vector containing six values ϕ i ¼ ½θ x θ y θ z x y z, where x is the x-coordinate and θ x is the angle of rotation around this axis.Analogous definitions are made for the y-axis and the z-axis.Here, S i ðϕ i Þ ¼ fs i;j g is the transformed skeleton (rotation and translation) and stands for the points at pose ϕ i E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 1 0 8 where R x ðθ x Þ is the rotation of angle θ x around the x axis.The pose estimation is based on the use of two successive views.The pseudocode is given in Algorithm 2. First, corners are detected in the region b i of the current frame F i using the FAST corner detector 26 with a nonmaximum suppression to obtain a better distribution of the corners C i ¼ fc i;j g over the MUAV, as illustrated in Fig. 4(a).The detected corners are described using the BRIEF 35 algorithm and matched with the corners from the previous frame, C i−1 , to obtain the matches M i;i−1 .At this point, the set of matches M i;i−1 contains corners located on the drone but also mismatches, corners on the background and corners on the propellers.In the next step, the unwanted corners (mismatch, background, and propellers) are filtered out by estimating the essential matrix 36 [see Fig. 4(b)].Hereafter M i;i−1 , C i−1 , and C i correspond to the filtered data set.
To estimate the pose of the drone, we need to associate the detected corners with points on the skeleton.A pinhole camera model is used to project the skeleton in the image plane.The camera is placed in the origin of the coordinate system with the standard axes, i.e., the depth of the scene is along the positive z-axis and the y-axis points to the ground.Using homogeneous coordinates, the camera projects the skeleton S i to S i;proj according to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 7 5 2 where P is the 4 × 4 homogeneous projection matrix. 37urther, in Fig. 4(c), we couple the detected corners and the projected skeleton points that are close together, dðc i−1 ; s i−1;proj Þ ≤ ϵ, c i−1 ∈ C i−1 and s i−1;proj ∈ S i−1;proj .If a skeleton point satisfies the condition with several corners, we keep the closest couple.From this we obtain subsets of C i−1 and S i−1;proj called C i−1;opti and S i−1;opti .It is then straightforward to obtain C i;opti , S i;opti , and couples ðc i;opti ; s opti Þ with c i;opti ∈ C i;opti and s opti ∈ S i;opti .
For each corner in C i;opti , a ray r passing through the camera center (the origin) and the corner is computed.The pose of the drone in image i, ϕ i , is then updated by minimizing the distance between each element of S i;opti and their corresponding rays (see Fig. 5).The minimization is performed using ϕ i−1 as a starting point E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 3 2 6 ; 5 5 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 5 0 4 where dðr j ; s j Þ is the orthogonal distance from the point s j to the line r j .Here, x c and y c are the corner coordinates in pixel.The minimization is performed by the Trust Region Reflective 38 algorithm, which is implemented in SciPy. 39ounds are used on the angular components of the pose to obtain meaningful results according to the MUAV specifications (see Table 2 in Sec. 6).The whole pose estimation algorithm is summarized as a pseudocode in Algorithm 2. The resulting pose vector ϕ i of the current frame and the state vector X i−1 of the previous frame are used to estimate the current state vector X i using the MUAV flight behavior model in Sec.6 and the Kalman filter in Sec. 5. Orientation (θ x , θ y , and θ z ) and velocity _ x are needed to estimate the current acceleration ẍi from the flight model [Eq.( 15)], taking into account the gravitation, thrust, and aerodynamic drag force.The determined acceleration vector and the ϕ i positions are then fed into the Kalman filter as measurement values (z i ) to update the MUAV state vector X i .Finally, the n-frame prediction is calculated as mentioned previously.
Figure 6 shows the results from the analysis of the same image sequence as used in Sec. 2. Owing to the missing reliable data for ground truth for the MUAV position, a reference for the 3-D flight path (-) is calculated from the state vectors (X i ).In Fig. 6(a), this track is compared to the five-frame prediction (gray bullets).A good agreement between flight path and prediction can be observed in area of constant or   linear movement (no rapid change of flight direction), while at turning points, where the flight direction changes rapidly, the algorithm still tends to overestimate the flight motion.In Fig. 6(b), a detailed analysis is presented.In the upper diagram, the amplitudes of the motion components, i.e., absolute velocity and absolute acceleration, are depicted.These values are directly derived from the algorithm.The velocity and acceleration are measured in standard units (m/s and m∕s 2 ).
In the lower diagram, the prediction error, that is L2-norm distance between predicted and ground-truth positions, is plotted for each frame.In most frames, this prediction error is < 20 cm.Only in a few sections of the sequence does this error exceed a value of 50 cm.These are the parts of the track where rapid changes of the flight direction occur.Further, significant prediction errors occur in the same frames as in the 2-D case [compare Fig. 2

(b)], only.
A more detailed discussion of the prediction error can be derived from Fig. 7. Here, the prediction error is plotted for each coordinate axis as the difference between the MUAV flight path (or track) and the predicted positions.Along the three coordinate axes, the prediction error is scattered around the zero position, giving evidence that no systematic error is present.In height (Δy) and width (Δx), the error is scattered with marginal standard deviation σ of AE1 and AE3 cm, respectively.Owing to a more dynamic flight behavior, the standard variation along the range axis (Δz) is σ ¼ 24.4 cm.Here, we want to point out that the retrieval of 3-D information from a single view in a 2-D intensity image is a challenging task, especially for range estimation.The statistical analysis of the whole sequence of 319 frames is summarized in Table 1.

Conclusions
We want to point out again that the 3-D information was carried out from pose estimation in 2-D intensity images.Further, in the analyzed sequences, the MUAV was observed under a very gradual observation angle.Here, a more slanted observation angle could help to increase the accuracy of the spatial position estimation.We successfully demonstrated the capability to extract 3-D information from 2-D intensity images by applying pose estimation with a simple skeleton model.From this pose information (i.e., position and orientation) and a simplified flight behavior model (with a floating body approximation), we could calculate the state vector that describes the current MUAV flight operation.Further, using a motion model that includes velocity and acceleration, we could predict the position of the MUAV in the near future.
The overall performance gives evidence that our approach is capable of resulting in a good first estimation of the future MUAV flight behavior.Unfortunately, in dynamic operation that includes rapid changes in velocity and direction, the underlying model tends to an overshooting effect.In future investigation, this error could be diminished by the use of a higher order of temporal derivative.
In our experiments, we used a single view from a single camera to retrieve a 3-D pose estimation by mainly analyzing the position of the four rotors.Further, the experiments were carried out at a very narrow viewing angle close to the MUAV principal plane.In these conditions, sometimes the line of sight on one or two rotors is blocked by the MUAV main body.Then, the pose estimation uses only the remaining visible rotors (three visible rotors are optimal).This situation is sufficient for a reliable pose estimation and has only a slight impact on the overall estimation of the 3-D position.In the later application in outdoor scenarios, we assume to have a steeper viewing angle, which should allow for an even better pose estimation.
To date, we have investigated our approach only in laboratory conditions as a proof-of-principle study.But, image data recorded in outdoor scenarios under operational conditions will contain a lot of additional disruptive factors such as signal-to-noise ratio, resolution, lighting conditions, and background.Further, MUAV will fly at larger ranges.Here, an active imaging device with a very narrow field of view could counter these challenges and could be used to gate out the background, increase resolution and signal-to-noise ratio, and be independent from lighting conditions.Our team has access to a shortwave-infrared laser gated-viewing system as presented by Christnacher et al. 10 In the near future, we will perform tests to investigate the pose estimation approach in operational conditions.Lighting conditions and observing MUAV flight behavior at a larger range can reduce the ranging accuracy of the proposed pose estimation approach due to an expected smaller scaling effect of the MUAV with the range at narrower observation angles.This effect could be compensated for by ranging capabilities of laser imaging devices.Notwithstanding this, the determination of the tilt angles (θ x and θ y ) and, thus, the application of the flight model should not be affected as long as the target is imaged with a high resolution.
Further, in the presented work, we have tested our approach on a single known MUAV, but we are convinced that our approach has a general nature and it could be easily adapted to other multicopter MUAV.Finally, the tracking and prediction of the flight behavior of unknown MUAV is possible if a classification of the target is done previously.As long as we have a proper skeleton model, this approach should be able to track and predict flight behavior of any multicopter MUAV.For the prediction of other MUAV types, such as fixed-wing or vertical takeoff and landing, the flight model has to be adapted carefully and in detail.
5 Appendix A: Brief Description of the Kalman Filter The Kalman filter is an iterative forward recursion process that consists of two steps: the "time update" and the "measurement update." 22Within the "time update," for the i'th measurement, the state vector X − i and the error covariance P − i are projected from previous values using the linear stochastic differential equation E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 5 6 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 5 2 2 Here, Eq. ( 10) reflects the process model as a system of motion equations with state-transition matrix A, the previous state vector X i−1 , and the control-input model Bu k (in our case, Bu k ¼ 0).Analogously, the projected error covariance matrix P − i is calculated in Eq. ( 11) using the error covariance P i and the process noise covariance matrix Q.
Then in "measurement update," the state vector X i is updated by weighting the measurement z i and the projected state vector X − i with the Kalman gain K i .Finally, the error covariance P i is updated E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 3 2 6 ; 4 0 2 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 3 5 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 3 2 6 ; 3 3 8 Here, the matrices H and R describe the measurement (or observation) process and covariance.
6 Appendix B: The Microunmanned Aerial Vehicles Flight Behavior Model The MUAV flight behavior can be described by the motion in Eq. (15).Here, m is the mass of the MUAV and mẍ is the resulting force that accelerates the MUAV in any direction.This force is the sum of the gravitational force G ¼ ½ 0 −mg 0 T , the total thrust T of all propellers, and the drag force D. Owing to the lack of information about the single motor thrusts, rotational motion is neglected in this model E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 1 8 3 The overall thrust T is assumed to be perpendicular to the MUAV main plane, that is, the plane in which all motors lie.Here, T is oriented in the direction of the normal vector n ¼ ½n x ; n y ; n z T .We define the thrust amplitude and the thrust vector as T ¼ jTjn ¼ ½T x ; T y ; T z T .The thrust amplitude has positive values jTj ∈ ½0; T max .Further, in a first approximation, we assume a floating body that means the vertical thrust (T y ) compensates for the gravitational force E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; From this assumption we can determine the forces induced by the thrust in the other direction E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 6 3 ; 6 7 Finally, we use a simplified model based on Stokes's law to determine the drag force x with k d as an aerodynamic constant.This force describes the aerodynamic drag on the MUAV that is proportional to the velocity.The aerodynamic constant k d can be determined experimentally. 40,41n this paper, we assume an isotropic aerodynamic coefficient of k d ¼ 0.55 kg s and a maximum thrust of T max ¼ 15.3 kg m s 2 , which can be calculated from Eq. ( 15) and the technical specification of the investigated MUAV. 34Details of this calculation are summarized in Sec. 7.
From Eqs. ( 15)-( 18), we can derive a calculus to determine the accelerations in the x-, y-and z-directions from the MUAV tilt angles (θ x and θ z ), which is analogous to the 3-D pose of the MUAV E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 6 3 ; 4 4 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 6 3 ; 3 8 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 6 3 ; 3 5 9 Further, stable flight operation has to comply with the condition that the thrust amplitude jTj is limited by the maximum thrust T max by jTj < T max .This condition has a direct impact on the expected values of the normal vector [see Eq. ( 22)] and is analogous to the definition of a maximum tilt angle θ max E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 6 3 ; 2 7 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 6 3 ; 2 1 0 the MUAV motion can be described as motion in the x-, yand z-directions by the equation system in Eqs. ( 25)-( 27) in analogy to Eq. ( 15).Now, θ i is the inclination of the MUAV or tilt angle in the i-direction and tan θ i is equal to n i n y E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 3 2 6 ; 5 2 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 3 2 6 ; 4 8 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 3 2 6 ; 4 6 1 In the following, the aerodynamic parameters T max and k d will be derived from this equation system using manufacturer specifications for weight, maximum horizontal and vertical velocity, and maximum tilt angle of the investigated MUAV. 34The values are summarized in Table 2.
As illustrated in Fig. 8, we can distinguish two borderline flight operations: maximum horizontal velocity and maximum vertical velocity.

Maximum Horizontal Velocity
First, we assume a flight operation with maximal horizontal velocity in the x-direction.In this case, the MUAV will be tilted to maximum inclination toward the flight direction, which is given by the maximum tilt angle θ ¼ θ max .Further, in the horizontal direction (x), the velocity _ x is equal to v h;max , the maximum horizontal velocity.
Further, due to borderline flight operation, the velocities in the other directions and all accelerations are zero, The thrust is at maximum level, T ¼ T max , but, due to the floating body assumption, its vertical component T y has to compensate for the gravitational force, that is T y ¼ mg.Thus, the maximal thrust T max can be calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 6 3 ; 6 9 6 T max ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi T 2 x þ T 2 y ; q E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 6 3 ; 6 4 5 ¼mg E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 ; 6 3 ; 6 1 6 T max ¼ 15;3 kg m s 2 : (30) Then, the aerodynamic constant k d can be calculated from Eq. ( 25) as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 1 ; 6 3 ; 5 7 3 k d ¼ mg tan θ max v h;max ¼ 0.55 kg s : Although we have already determined the values for both T max and k d , these results can be validated by the second case, that is, maximum vertical velocity.

Maximum Vertical Velocity
Here, we can assume a flight operation with maximum vertical velocity ẏ ¼ v v;max .In this case, the inclination angle is zero, θ ¼ 0 deg, and the maximum thrust pushes the MUAV upward in the y-direction, T y ¼ T max .The velocities in the horizontal directions are zero, _ x ¼ _ z ¼ 0, as well as any accelerations, ẍ ¼ ÿ ¼ z ¼ 0.Then, from Eq. ( 26) and with the result for either k d or T max of the previous case, we can again derive E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 2 ; 6 3 ; Both results agree with the previously found values.

Algorithm 1
Fig.1CMT of the key points in a MUAV video sequence.

3
Prediction of Microunmanned Aerial Vehicles Position in Three-Dimensional Space In three-dimensional (3-D) space, the flight path prediction has to take into account the position and motion in all spatial dimensions to reflect the complex MUAV flight behavior.Here, we want to introduce a pose estimation algorithm to derive the 3-D position and orientation of the MUAV from the 2-D intensity images.Further, we use the obtained data in an aerodynamic flight behavior model to calculate the acceleration as an additional measurement value.Then, both the determined 3-D position and the acceleration values are fed into the aforementioned Kalman filter and prediction algorithm to return the probable future MUAV position in the 3-D space.

Fig. 2
Fig. 2 Prediction of MUAV flight behavior in the 2-D image plane using the CMT tracker and the first and second derivatives ( _ x and ẍ ): (a) ground-truth CMT track (-) and five-frame track prediction using _ x -only (Δ) and _x -ẍ (□) algorithms, (b) the amount of the nj _ x j (⋯) and n 2 j ẍ j (-) components and the prediction errors that are the L2-norm distances from the original track.

Fig. 3
Fig. 3 Representation of a drone as (a) a computer graphics model and (b) a skeleton, as used in the pose estimation of the quadcopter.

Algorithm 2 1 : 1 3:
Pseudocode for the pose estimation of the MUAV from two successive views.Input: Current frame F i , previous pose ϕ i−1 ∈ R 6 , previous detected corners C i−1 , b i from the CMT, threshold ε, skeleton S Output: ϕ i , C i Detect corners C i within b i using the FAST algorithm 2: Match C i with C i−1 using BRIEF descriptors and Hamming distance → M i;i−Filter out propellers and background motion from M i;i−1 by estimating the essential matrix 4: Use

Fig. 4
Fig. 4 Link between the detected corners on the image and the drone skeleton model.(a) The detected corners using FAST, (b) the corners after the estimation of the essential matrix and the projected skeleton from the previous frame, and (c) skeleton points and corners used for the pose estimation of the drone.

Fig. 6
Fig. 6 Results of the 3-D MUAV tracking and n-frame prediction algorithm applied to a video sequence recorded under laboratory conditions: (a) 3-D tracks from the Kalman filter (-) and the n-frame prediction (gray bullets); (b) amplitude of the components (first and second derivatives) and the prediction mismatch (L2-norm).

Fig. 5
Fig. 5 Principle for aligning the drone skeleton on the rays issued from the detected corners.Aligning the skeleton S starting at the pose ϕ i−1 on rays issued from C i;opti shown in Fig. 4.

Fig. 7
Fig. 7 Illustration of the prediction error as a 3-D scatter plot.The difference between track and prediction is plotted for each coordinate axis: height, width, and range.

) 7
Appendix C: Estimation of Aerodynamic Coefficient k d and Maximum Thrust T max With the definitions of the gravitational force G, the thrust T, and the drag force D as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 6 3 ; 1 2 3 G ¼

Fig. 8
Fig.8The two different extreme flight operations, case 1: maximal horizontal velocity and case 2: maximal vertical velocity, are distinguished.
T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 3 ;6 3 ; 3 2 4 T max ¼ mg þ k d v v;max ¼ 15.3 kg m s 2 :

Table 1
Summary of statistics on the prediction error in a sequence of 319 frames.

Table 2
34ysical flight parameter T max and k d for the investigated MUAV calculated from the manufacturer specifications.34