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

UnmannedAerial 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.


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 2 International Journal of Aerospace Engineering 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.

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 selfaware 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 ∈ [1, 2, . . . , ] to indicate its position along the branch. Thus, the coordinates of the antennas can be written as where is the Kronecker delta, that is, 1 for = and 0 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 as where 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: 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 ∼ (0, 2 ) is a random uniform variable between 0 and 2 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.

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 aŝ( In [8] a technique to estimatê( ) 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 International Journal of Aerospace Engineering 3 (2).̂( ) are estimated as ( ) , where ⋅ stands for the phase operator, where where an average over /2 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 arctan(⋅, ⋅) is the arctangent function extended to the [0, 2 ) 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.

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 is while 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 { ( ) } =1,..., and a set of weights { ( ) , } =1,..., . 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] { ( ) −1,RES } using the regression constant m −1 : where ( ) −1,RES is the th particle of − 1 step cloud after the resampling process [11]. ∼ N ( , = 0) is a Gaussian random variable that represents the fundamental statistical scattering part of the motion model. The parameter m −1 is the slope of the linear regression of previous values oft hat 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 2 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 2 offsets of | ( ) −̃| due to the periodicity of , the minimum absolute value of the distance of the two angles ( ) and̃in 2 -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 ( ) , ,RESH is the weight after the reshaping process. Finally, after the estimation stage, the particles cloud is resampled, generating the set ( ) ,RES , to eliminate spurious particles and avoid fragmentation.

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 Ω ≈ 6 ∝ s . Likewise the PC-G can be estimated as ≈ 6 /( −1) ∝ . 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 estimatioñintroduced in (5) and the resampled particle swarm ( ) ,RES as follows: where the operators avg and std 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 > 0 is the proportional feedback parameter of the PSAS. The value of , depending on the two feedback parameters + and − and the threshold value th , 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 [− th , th ] convergence interval, that is, when − > th or + < − th . 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 [− th , th ] convergence interval, that is, when + > th or − < − th . 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 th are related to the dynamic behaviour of the convergence point that is related to the UAV motion model. Optimal values for and th were numerically calculated for the trajectories generated by CARATE in this paper.

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.
International Journal of Aerospace Engineering 5 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 as The angle , depicted in Figure 1, is the local azimuthal angle at time of the actual trajectory coordinates [ , , ] seen from the previous path point [ −1 , −1 , −1 ]. The local elevation angle can be defined analogously. V is the absolute speed of the UAV. At each step, the trajectory increments generated by , , and V are related to the "seed" components seed, , seed, , and V seed, and to the random components , , and V as follows: = seed, + , = seed, + , 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 V > 0, the UAV will have an accelerating trend and for > 0 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 by (17) In 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 ang and sp , 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 V of the random components , , and V introduced in (15).
The CARATE algorithm, as visible in the sample path in Figure 2, generates a 3D trajectory inside a limited flight area, formed by the horizontal bound region and the vertical bound region .
is defined as {| | < ∩| | < } while is defined as { > min ∩ < max }. A specific routine of the algorithm has been developed to keep the trajectory inside the limited flight area, preserving it smooth and without sharp edges on region bounds. If or exceed the value of seed, introduced in (16) is forced for the subsequent bend steps to follow a set of equispaced values between seed, and seed,END . The value of seed,END is appropriately defined depending on the Cartesian quadrant where the trajectory crosses the surface of . Analogously, if exceeds the value of seed, introduced in (16) is forced for the next bend steps to a set of equally spaced values between seed, and 6 International Journal of Aerospace Engineering seed,END . The value of seed,END is also defined depending on the Cartesian quadrant where the surface of has been crossed.

Cramér-Rao Lower Bound
In statistics the Cramér-Rao Lower Bound (CRLB) expresses the lowest value of the Root Mean Square Error (RMSE) of an optimal estimator given a certain set of data [13]. In this paper the estimations of the coordinates [ , ] are carried out from the information given by the measurements ( ) , ( ) , and ( ) . , ] as the vector of ideal measured values, namely, ( ) , in (4) but without the detrimental effects due to the impairments introduced in (2) such as PN, DS, and C-AWGN: where, as introduced in Section 2,̂( ) is the D-PoA for th 3A-ULA antenna branch and | | is the constant amplitude of the complex downconverted constellation. The vector of real measurements u = [ ( ) , ( ) , ( ) ], defined in (4), can be written as where the intermediate signal | | 2̂( ) + ( , ) is averaged through antenna pairs and̂= (2/ ) ∑ where the notation (⋅) * denotes the complex conjugate operator. ( , ) is the C-AWGN with zero mean and variance 2 . ( , ) is the PoA as defined in (3). The expectation of the differential noise ( , ) among the antennas of the same branch is = 0, while its variance 2 = 2| | 2 2 + 4 .
To calculate the CRLB it is necessary to evaluate the loglikelihood ( , ) = ln [u | b ] of the acquired data given the parameters to be estimated. ln(⋅) represents the natural logarithm. In this context, for the calculation of CRLB, the measurements vector u is considered affected only by C-AWGN noise without the influence of the other impairments introduced in (2). This assumption yields to the likelihood Gaussianity and is verified via numerical fitting similarly to [14], so that where the time index is omitted until the end of the CRLB derivation for notational simplicity. The CRLB of the variable , defined as th element of the unknowns vector = [ , ], is where ( , ) is the Fischer information matrix: For the log-likelihood function described in (21), we obtain where = [cos sin , sin sin , cos ] is the unit direction vector and norm = / 0 is the normalized interantenna distance (Figure 4). The notation (⋅) T denotes the transpose operator. Finally, the CRLB expressions for and are The RMSE of every estimator of the unknown from the measurements u is lower bounded by √CRLB . It is important to state again that this C-AWGN-only CRLB calculated in (25) takes into account C-AWGN but not the other previously introduced impairments, in particular DS and PN. Given the previous assumptions, the AWGN-only CRLB is lower or equal to the fully impaired CRLB. Another relevant point is that the computed CRLB does not take into account the motion model and the trajectory history differently from the PF algorithm. The CRLB is simply computed from the 3 observations ( ) introduced in (4) at time , as the BM does. Instead, the proposed PF algorithm exploits the motion model jointly with PSAS, using information coming from the trajectory history. For this reason, it may occur that the PF estimation offers better performance than the CRLB.

Performance Analysis
The performance of the PF algorithm is now assessed. A comparison with the BM estimation introduced in Section 3 and the CRLB calculated in Section 7 is also made. Different UAV trajectories are generated with the CARATE to assess performance. Each UAV trajectory is time samples long. In order to have a common performance evaluation metric, in are the RMSEs calculated for each trajectory traj . traj is the number of generated trajectories. Analogously, with the purpose to allow a performance comparison with the CRLB, we define the average trajectory CRLB (AT-CRLB), for ∈ [ , ] and where The default parameters for the performance evaluation are (a) angles = 7 ∘ , = 7 ∘ , and ang = 20; (b) speeds V min = 0.5 km/h, V max = 10 km/h, and sp = 2. The array parameters are = 6 and norm = 1/3. The PF algorithm basic settings are = 20 and = 20. PSAS parameters = 25% and th = 10 ∘ were numerically optimized for the CARATE algorithm. In the signal model introduced in (2), the PN is considered white with standard deviation PN = 1 ∘ , the PO is considered uniform in [0, 2 ) for each array element, and the signal-to-noise ratio is SNR = [| | 2 / 2 ] dB =10 dB. The number of simulated trajectories traj and their length have been properly dimensioned to guarantee the statistical confidence of the results.
In Figure 2 an example of 3D-trajectory generated for performance evaluation is shown. The 3A-ULA is positioned in the center of the axis.

Performance with Static and Dynamic
. We now consider the performance as a function of the signal-to-noise ratio SNR = | | 2 / 2 . In Figure 3, the AT-RMSE of the azimuthal coordinate and the elevation coordinate for both the PF algorithm and the BM are reported. Furthermore the impact on performance of the PSAS algorithm introduced in Section 5 is analysed. It is possible to observe how, due to the geometrical asymmetry, the performance in terms of AT-RMSE is significantly different between and . Greater SNR values lead to lower AT-RMSE for both PF and BM estimations. It is important to state that, for a static value of , PF performs better than BM in a broad range of SNR values. For example, at 2 dB of SNR, it outperforms the BM algorithm by about 17 ∘ for and 10 ∘ for . The analysis shows that the introduction of the PSAS algorithm leads to a further improvement of performance. For example, at 2 dB of SNR the PSAS leads to a performance improvement of 5 ∘ for and 3 ∘ for with respect to the use of a static . It is visible in Figure 3 how the introduction of the PF algorithm allows overcoming the AT-CRLB. This is possible thanks to the a priori knowledge of the UAV path extracted by the motion model from the previous steps of the trajectory. For higher SNR the performance improvement of PF with respect to the BM algorithm is less prominent than for lower SNR.

Performance for Different Interantenna Distances norm .
Now, the behaviour of the AT-RMSE for the BM and PF algorithms is analysed varying the normalized interantenna distance norm . To unambiguously extract the phase from ( ) the interelement distance must be less than half of the impinging signal wavelength [15]: norm ≤ 0.5. Furthermore in (25) it is visible that for higher values of norm the CRLBs increase. This leads to norm = 0.5 as the optimal value for the estimation. However, in lower SNR scenarios and for norm values near 0.5 the CRLB is found to be overly optimistic. In fact, high C-AWGN and the other impairments in (2) severely impair performance corrupting the value of ( ) in (4) and making it exceeding 180 ∘ . The PF algorithm, thanks to its hypothesised a priori knowledge of the UAV path given by the motion model, attenuates this impairing effect. In Figure 3, we can see that the proposed PF algorithm is less dependent from norm = 0 / 0 than the BM. Thus, the same 3A-ULA is usable for a wider set of carrier frequency 0 with respect to the BM that exhibits its optimal performance for a narrow interval of norm around 0.4. For this reason, the PF algorithm can offer good performance also in a multiple UAVs scenario where multiplexing is implemented in a frequency division fashion.

Conclusions
We have discussed the application of an appropriately designed PF algorithm for self-localization (SL-S) and ground localization (GL-S) of UAVs using a 3A-ULA of antennas, showing an overall increase of performance with respect to the BM. A novel algorithm to manage the amplitude of the particle swarm, namely, Particles Swarm Adaptive Scattering (PSAS), has been developed and tested, showing a further increase of precision. A complete, fully adjustable, and effective 3D UAVs trajectory emulator, namely, CARATE, has been proposed and used to assess performance. Strong impairing effects like Doppler spread, phase offset, and phase noise have been considered in the performance evaluation. The effects on the proposed localization algorithm of the PF model parameters as well as the SNR and the 3A-ULA characteristics have been studied. Numerical results show that the proposed PF algorithm is able of dynamically track the UAVs angular position better that the BM. A critical point of this approach is that, similarly to PSAS behaviour, a dynamic calibration procedure has to be implemented to optimize the PF parameters. Future work will also investigate