International Journal of Aerospace Engineering

Volume 2016 (2016), Article ID 7630950, 9 pages

http://dx.doi.org/10.1155/2016/7630950

## Azimuth and Elevation Dynamic Tracking of UAVs via 3-Axial ULA and Particle Filtering

^{1}WiPLi Lab, University of Udine, Via delle Scienze 206, 33100 Udine, Italy^{2}EcoSys Lab, Alpen-Adria-Universität Klagenfurt, Universitätsstraße 65-67, 9020 Klagenfurt, Austria

Received 1 April 2016; Revised 23 May 2016; Accepted 29 May 2016

Academic Editor: Yanbo Zhu

Copyright © 2016 Andrea Papaiz and Andrea M. Tonello. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Unmanned Aerial Vehicles (UAVs) localization has become crucial in recent years, mainly for navigation or self-positioning and for UAV based security monitoring and surveillance. In this paper, azimuth and elevation radio positioning of UAVs are considered. The localization is based on multiple differential phase-of-arrival measures exploiting a 3-Axial Uniform Linear Array of antennas. An ad hoc particle filtering algorithm is applied to improve the positioning performance using a dynamic motion model. A novel adaptive algorithm, namely, Particles Swarm Adaptive Scattering (PSAS), is proposed to increment the algorithm stability and precision. To assess performance a Confined Area Random Aerial Trajectory Emulator (CARATE) algorithm has been developed to generate actual paths of flying UAVs. The algorithm performance is compared with the baseline method and with the average trajectory Cramér Rao lower bound to show the effectiveness of the proposed algorithm.

#### 1. Introduction

Unmanned Aerial Vehicles (UAVs) are attracting considerable attention since they can be used for a number of consumer, industrial, and military applications ranging, for instance, from sport video making to environmental monitoring and parcel delivery [1, 2]. A key component is the technology that allows to locate and navigate the UAVs. Presently, global navigation satellite systems (GNSS) and inertial sensors are used to provide information on position, speed, and direction of movement. Reliable localization is very important also in view of new regulations that aim at better controlling the use and the status of the UAVs for higher safety and security [3]. Furthermore, recently, UAVs technology has started to be under the spotlight of police audits because of possible security threats [4].

In this paper, a radio localization approach is considered and it is based on azimuth and elevation positioning using a transmitting source as a reference. Azimuth and elevation are determined by processing with particle filtering (PF) the signals that impinge on a 3-Axial Uniform Linear Array (3A-ULA). The 3A-ULA can be mounted either on a ground base station or on the UAVs. In the first case, namely, Ground-Localization Scenario (GL-S), the base station passively eavesdrops the signals emitted by the UAVs to determine their angular coordinates. In the latter case, namely, Self-Localization Scenario (SL-S), the ground node acts as a radio anchor allowing UAVs self-localization. The system can be used in a standalone way or to complement existing GNSS/inertial or inertial sensors employing data fusion techniques [5].

In Section 2, the analytical model for the signals received by the 3A-ULA is described. The baseline method (BM) coordinates estimates obtained from the 3A-ULA in Section 3 are iteratively processed with an ad hoc particle filtering algorithm (PF) described in Section 4. The PF technique uses a novel dynamic particles swarm management routine, namely, Particles Swarm Adaptive Scattering (PSAS), introduced in Section 5 to improve both convergence and precision performance.

In order to quantify the performance of the proposed tracking algorithm an emulator for the UAV behaviours is specifically developed in Section 6. This ad hoc numerical model, namely, Confined Area Random Aerial Trajectory Emulator (CARATE), randomly generates trajectories in the 3-dimensional space following a well defined set of rules in order to emulate real UAV behaviours.

Finally, in Section 8, numerical results are reported to assess performance as a function of the UAVs motion model, the PF algorithm parameters, and external detrimental impairments, such as phase noise (PN) and carrier frequency Doppler shift (DS). A comparison with the average trajectory Cramér Rao lower bound (AT-CRLB) for the angular estimation of the elevation and azimuth is also provided.

#### 2. Positioning System Description

A system that determines the position of flying UAVs through the analysis of a signal impinging on the 3 branches of a 3A-ULA is considered.

The GL-S is built as follows. UAVs send narrow band radio signals that are captured by a ground base station equipped with a 3A-ULA. The radio signals are properly frequency/time duplexed to allow multiple UAVs tracking. The ground node can then compute locally the UAVs azimuth and elevation coordinates [6]. By sharing the data collected from more ground nodes, it is possible to achieve a full 3D positioning for the UAVs trajectories.

The SL-S is analogously arranged. Each 3A-ULA is mounted on the UAVs so that the ground device takes the role of a radio beacon anchor [7], allowing the UAVs to be self-aware of their position. Using more anchors transmitting at different frequency or in a time division multiplexing fashion, the UAVs can accomplish a full 3D positioning. Clearly, the implementation of the SL-S is dependent on the UAV and array sizes.

The 3A-ULA is constituted, for each one of the 3 branches, by antennas. Each branch is labelled with to indicate along which one of the axes it is displaced. An antenna of the same branch is indexed with to indicate its position along the branch. Thus, the coordinates of the antennas can be written aswhere is the Kronecker delta, that is, for and otherwise, and is the constant distance between antennas.

The signal model considered in the following describes for simplicity the scenario with only one UAV; this assumption can be easily extended to a multi-UAVs application. Both GL-S and SL-S scenarios, due to symmetry, can be modelled in the same way. Hence, in general, the downconverted signal received by the th antenna of the 3A-ULA branch can be written aswhere indexes the time, sampled with period . The quantity is the phase of the demodulated signal impinging on the th sensor of the branch , namely, the phase of arrival (PoA), in the ideal case without impairing effects. represents the phase Doppler shift (DS) caused by the nonzero speed of the radio source with respect to each receiving array element. The parameter is the phase noise (PN) of the oscillators. It is considered white for this application [8]. The component is the phase offset (PO) between modulating and demodulating oscillators. The phase components of (2) are described in more detail in the following:where is the carrier frequency of the received modulated signal and where is the time of arrival (ToA) of the plane wave from the source to the sensor. is the DS coefficient due to the movement of the UAV. The parameter is the UAV speed in the direction of the center of the 3A-ULA and is the speed of light. Finally, is circular additive white Gaussian noise (C-AWGN) with standard deviation . The receiving architecture comprises downconverting branches for each axis of the 3A-ULA all fed by the same oscillator. For this reason PN and PO components are identical for all antennas.

Considering the scenario where the transmitted signals carry also information, the component of (2) is representative of the unknown informative part of the signal, specifically the modulated complex data symbols. For the application considered in this paper is assumed to be an equiprobable symbol of an unknown order PSK modulation scheme. For this reason the symbol is modelled as complex uniform ring shaped noise (RSN) where is a random uniform variable between and and is the constant amplitude of the modulated signal.

In this model, the presence of multipath propagation is not considered; a line of sight environment is assumed. This can hold true when the ground and UAV antennas have a wide vertical lobe, respectively, directed upward and downward.

#### 3. Azimuth and Elevation Estimation with 3A-ULA

A radio source moving in the 3D Cartesian space is considered. Its angular spherical coordinates and are determined using the same coordinate system of the 3A-ULA defined in (1). A trajectory example is visible in Figure 2. The positioning technique used to locate the radio source is based on the estimation of the phase difference of the received signals between different antenna couples belonging to the same 3A-ULA branch , namely, the differential phase of arrival (D-PoA) . The D-PoA is defined as .

In [8] a technique to estimate in a scenario with impaired signals similar to (2) is described. Applying this technique for the model considered in this paper, it allows exploiting the D-PoAs from the signal model introduced in (2). are estimated as , where stands for the phase operator, wherewhere an average over distinct antenna pairs is performed to mitigate noise and compensate other zero mean impairments. Using , the angular spherical coordinates and are estimated, respectively, with and as follows:where is the arctangent function extended to the domain.

The estimations computed in (5) will be referred to as baseline method (BM) estimations to differ them from the particle filter (PF) estimations that will be introduced in Section 4.

#### 4. PF Algorithm Estimation

In this section, the application of the PF algorithms [9] to the localization problem introduced in Section 3 is briefly described. PF algorithm techniques have been chosen because of their simple, scalable, and flexible numerical implementation with respect to other methods like Extended Kalman Filter [10] that approach the problem linearising the system and approximating the probability density functions as Gaussians. PF is specifically designed to fit complex system models as well as nonlinear and non-Gaussian configurations. The target of PF methods is to estimate the dynamic evolution of the hidden states in a system that are the coordinates and of the radio source, using the sequential discrete noisy measures introduced in (4). In PF, it is necessary to hypothesize a model for the observations and for the states temporal evolution [9]. The evolution model for the hidden states and iswhile the model for the observations related to the hidden states introduced in (4) is

The relations (6) show how the current hidden state evolves thanks to the previous states and the random uncertain variables and through the motion model function . Relation (7) shows how the measurements depend only on the current state and on a random uncertain variable through the measure model function .

For notational simplicity, both and will be denoted with the generic angle . When needed this generalization will be removed.

PF has been chosen because of its simplicity and versatility. PF algorithms are a numeric implementation of Bayesian estimation [10], a statistical approach to iteratively increase the knowledge of the states to be estimated by narrowing their pdf conditioned by the sequential acquisition of new measures during time. The pdf of , defined as , is approximated in a discrete way [11] with a set of particles and a set of weights . Particles and weights are generated from the partial knowledge given from the measure model and the state evolution model . This knowledge is more and more incremented by the continuous time acquisition of the measures . Clearly, the UAV moves so that its coordinates change over time .

As shown in (6), in the PF algorithm an assumption about the model that controls the evolution of the hidden state of the system is made. For this purpose the motion model used is the Drift Motion Model (DMM) proposed in [12]. The DMM generates the new particles cloud at step shifting the previous step resampled cloud [9] using the regression constant :where is the th particle of step cloud after the resampling process [11]. is a Gaussian random variable that represents the fundamental statistical scattering part of the motion model. The parameter is the slope of the linear regression of previous values of that are the PF estimates of . It is important to emphasize that in this regression model, due to the periodic nature of , it is necessary to use the previously saved unwrapped values of . This is necessary to correctly track the azimuthal coordinate during angular shifts from to going from the third to the second Cartesian quadrant.

Particles weights are calculated by the PF algorithm using in each cycle the temporary estimates obtained from the real acquired measures in (5). This temporary estimations are compared with the cloud of generated particles through an exponential weighting function to obtain the set of weights :where, to correct unnecessary offsets of due to the periodicity of , the minimum absolute value of the distance of the two angles and in -modulo space is computed. A reshaping technique [12] is then applied, before estimation, to regularize the particle cloud weights.

The PF estimation of the hidden state is then computed through a weighted average of the reshaped particles swarm:where is the weight after the reshaping process. Finally, after the estimation stage, the particles cloud is resampled, generating the set , to eliminate spurious particles and avoid fragmentation.

#### 5. Novel Particles Swarm Adaptive Scattering Algorithm (PSAS)

The parameter is the standard deviation of the random variable introduced in (8). In the previous section , it was considered constant over time. However, depending on the situation of the tracking algorithm, it will be shown that it is better to have a dynamic . It is possible to define two parameters that describe the particle cloud: , namely, the Particles Cloud Width (PC-W) that describes how much wide the set of particles is and, , the Particles Cloud Granularity (PC-G) that indicates the interparticle distance:The value of influences the convergence speed of the algorithm. In fact, if the particles cloud is positioned over the value to be estimated the convergence is fast because the the exponential weighting function introduced in (9) gives high weights to the nearest particles and negligible weights to the further ones. Instead, if the particles cloud is located far from the value to be estimated, convergence is slower because the exponential weighting function algorithm gives almost the same weights to all particles. Thus, it will need more iterations to move the particles swarm, cycle by cycle, near the convergence point. For these reasons a bigger value of will lead to wider clouds, encouraging the convergence. On the other hand, smaller values of will make the algorithm slower and more unstable to fast state changes.

The value of influences the precision of convergence. In fact, hypothesising that it has reached the convergence point, the estimation precision is upper bounded by the average interparticle distance, that is, the PC-G.

Considering the spreading variable introduced in (8) as a Gaussian variable with zero mean and standard deviation , it is possible to estimate the PC-W using the Tchebysheff’s inequality as . Likewise the PC-G can be estimated as . These relations show that is directly proportional to PC-W and PC-G. Because of the previous observation it is possible to conclude that, given a swarm of particles, during the convergence phase a bigger value of would be more appropriate because it would lead to faster results. On the other hand, after the convergence phase, smaller values of would be more appropriate to increase the estimation precision.

For these reasons in this paper a novel algorithm that manages the dynamic evolution of across the iterations, namely, Particle Swarm Adaptive Scattering (PSAS), is proposed.

The distance metric is defined as a baseline error cloud for the particles and is used by the PSAS to select the convergence status of the PF algorithm. The error cloud is calculated between the baseline estimation introduced in (5) and the resampled particle swarm as follows:where the operators and denote, respectively, the numeric operations of average and standard deviation and the absolute value operates in the periodic domain. The parameters and represent, respectively, the lower and higher bounds of the error cloud and are used as thresholds by the PSAS algorithm to trigger the increment or the decrement of as follows:where is the proportional feedback parameter of the PSAS. The value of , depending on the two feedback parameters and and the threshold value , is set positive to iterative increment and negative to decrement and is set to zero to keep it constant.

The value of is incremented with respect to its previous value when the condition of nonconvergence is detected. This condition holds true when the error cloud is reasonably out of convergence interval, that is, when or . Then, the PC-W is incremented to encourage convergence. If the previous condition is not satisfied, the value of is decremented w.r.t. its previous value when the convergence state is detected. This condition holds true when the error cloud is on the edges of convergence interval, that is, when or . Then, the PC-G is reduced to increment the precision. For the remaining values of and the condition of tight convergence is reached. In this state the value of is kept the same w.r.t. its previous value. This PSAS case keeps the cloud at the same size to avoid instability conditions due to rapid changes of the convergence point. The values of the feedback parameter and the threshold parameter are related to the dynamic behaviour of the convergence point that is related to the UAV motion model. Optimal values for and were numerically calculated for the trajectories generated by CARATE in this paper.

#### 6. UAVs Trajectory Emulator

To assess performance of BM and PF algorithms estimations and to evaluate the effects of different parameters on the estimation process, it is necessary to emulate UAVs trajectories using a proper coordinates evolution algorithm. Herein, a possible trajectory generation method is proposed.

The proposed algorithm, named Confined Area Random Aerial Trajectory Emulator (CARATE), generates iteratively a 3D path obtained from a variable length previous history of the trajectory and a tunable set of random variables. CARATE is specifically designed to emulate UAVs trajectories inside a limited flight area, in the GL-S, around the receiving array. Thus, the generated trajectories allow testing the tracking algorithms for a broad range of arrival angles. An UAV located far from the receiving array instead would be localized using a limited range of angles. The Cartesian position at time of the UAV evolves asThe angle , depicted in Figure 1, is the local azimuthal angle at time of the actual trajectory coordinates seen from the previous path point . The local elevation angle can be defined analogously. is the absolute speed of the UAV. At each step, the trajectory increments generated by , , and are related to the “seed” components , , and and to the random components , , and as follows:The random components are generated following, respectively, the normal distributions , , . The parameters of the random component distributions affect the UAV path behaviour. For example, for , the UAV will have an accelerating trend and for the UAV will tend to turn anticlockwise. The seed values introduced in (15) are the history parameters that express how much the evolution of the path is related to its past and are given bywhere is the arctangent function defined within and whereIn Figure 1, the relations among the CARATE components for the variable are represented. As shown in (16), the smaller the window width of the seed parameters and , the higher the influence of the path randomization with respect to the past history of the path itself. The same effect will occur for higher values of the standard deviation , , and of the random components , , and introduced in (15).