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

Authors: Andrea Papaiz; Andrea M. Tonello
Journal: International Journal of Aerospace Engineering
 (2016)

---

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


---

## Body

## 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 Section2, 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 Section6. 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 Section8, 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, byNa antennas. Each branch is labelled with a∈x,y,z to indicate along which one of the axes it is displaced. An antenna of the same branch a is indexed with i∈[1,2,…,Na] to indicate its position along the branch. Thus, the coordinates of the antennas can be written as(1)xa,i=i-1Dδxa,ya,i=i-1Dδya,za,i=i-1Dδza,where δij is the Kronecker delta, that is, 1 for i=j and 0 otherwise, and D 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 theith antenna of the 3A-ULA branch a can be written as(2)sna,i=Anejψna,i+γna,i+φn+ζ+wna,i,where n indexes the time, sampled with period T. The quantity ψn(a,i) is the phase of the demodulated signal impinging on the ith sensor of the branch a, namely, the phase of arrival (PoA), in the ideal case without impairing effects. γn(a,i) represents the phase Doppler shift (DS) caused by the nonzero speed of the radio source with respect to each receiving array element. The parameter φn 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:(3)ψna,i=2πf0τna,i,γna,i=2πf0CDSnT-τna,iτna,i,where f0 is the carrier frequency of the received modulated signal and where τn(a,i) is the time of arrival (ToA) of the plane wave from the source to the sensor. CDS(t)=vr(t)/c0 is the DS coefficient due to the movement of the UAV. The parameter vr(t) is the UAV speed in the direction of the center of the 3A-ULA and c0 is the speed of light. Finally, wn(a,i) is circular additive white Gaussian noise (C-AWGN) with standard deviation σw. The receiving architecture comprises Na 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 componentAn 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 An 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) An=|A|ejξn where ξn~U(0,2π) is a random uniform variable between 0 and 2π and |A| 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θn and ϕn 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 a, namely, the differential phase of arrival (D-PoA) ψ^n(a). The D-PoA is defined as ψ^n(a)=ψn(a,2i-1)-ψn(a,2i)∀i∈[1,2,…,Na/2].In [8] a technique to estimate ψ^n(a) 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). ψ^n(a) are estimated as un(a), where · stands for the phase operator, where(4)una=2Na∑i=1Na/2sna,2i-1sna,2i∗,where an average over Na/2 distinct antenna pairs is performed to mitigate noise and compensate other zero mean impairments. Using un(a), the angular spherical coordinates θn and ϕn are estimated, respectively, with θ~n and ϕ~n as follows:(5)θ~n=arctan⁡unx2+uny2,unz,ϕ~n=arctan⁡uny,unx,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.


## 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 θn and ϕn of the radio source, using the sequential discrete noisy measures un(a) 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 θn and ϕn is(6)θn=gθθn-Nmm÷n-1,pθ,n,ϕn=gϕϕn-Nmm÷n-1,pϕ,n,while the model for the observations related to the hidden states introduced in (4) is(7)una=hθn,ϕn,qn.The relations (6) show how the current hidden state [θn,ϕn] evolves thanks to the previous Nmm states and the random uncertain variables pθ,n and pϕ,n through the motion model function g. Relation (7) shows how the measurements un(a) depend only on the current state [θn,ϕn] and on a random uncertain variable qn through the measure model function h.For notational simplicity, bothθn and ϕn will be denoted with the generic angle αn. 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 αn, defined as fαnα, is approximated in a discrete way [11] with a set of M particles {αn(m)}m=1,…,M and a set of weights {wα,n(m)}m=1,…,M. Particles and weights are generated from the partial knowledge given from the measure model h and the state evolution model gα. This knowledge is more and more incremented by the continuous time acquisition of the measures un(a). Clearly, the UAV moves so that its coordinates change over time n.As shown in (6), in the PF algorithm an assumption about the model that controls the evolution of the hidden state of the system gα 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 {αn(m)} at step n shifting the previous step resampled cloud [9] {αn-1,RES(m)} using the regression constant mn-1:(8)αnm=αn-1,RESm+mn-1T+η,where αn-1,RES(m) is the mth particle of n-1 step cloud after the resampling process [11]. η~N(σs,ms=0) is a Gaussian random variable that represents the fundamental statistical scattering part of the motion model. The parameter mn-1 is the slope of the linear regression of Nmm previous values of α^n that are the PF estimates of αn. It is important to emphasize that in this regression model, due to the periodic nature of αn, it is necessary to use the previously saved Nmm unwrapped values of α^n. This is necessary to correctly track the azimuthal coordinate ϕn during 2π angular shifts from -π to +π going from the third to the second Cartesian quadrant.Particles weightswα,n(m) are calculated by the PF algorithm using in each cycle the temporary estimates α~n obtained from the real acquired measures in (5). This temporary estimations are compared with the cloud of generated particles αn(m) through an exponential weighting function to obtain the set of weights wα,n(m):(9)wα,nm=exp-αnm-α~n,where, to correct unnecessary 2π offsets of |αn(m)-α~n| due to the periodicity of αn, the minimum absolute value of the distance of the two angles αn(m) and α~n in 2π-modulo space is computed. A reshaping technique [12] is then applied, before estimation, to regularize the particle cloud weights.The PF estimationα^n of the hidden state αn is then computed through a weighted average of the reshaped particles swarm:(10)α^n=∑m=1Mwα,n,RESHmαnm,where wα,n,RESH(m) is the weight after the reshaping process. Finally, after the estimation stage, the particles cloud is resampled, generating the set αn,RES(m), to eliminate spurious particles and avoid fragmentation.


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

The parameterσs is the standard deviation of the random variable η introduced in (8). In the previous section σs, 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 σs,n. It is possible to define two parameters that describe the particle cloud: Ωn, namely, the Particles Cloud Width (PC-W) that describes how much wide the set of particles is and, εn, the Particles Cloud Granularity (PC-G) that indicates the interparticle distance:(11)Ωn=maxmαnm-minmαnm≤2π,εn=ΩnM-1.The value of Ωn 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 Ωn will lead to wider clouds, encouraging the convergence. On the other hand, smaller values of Ωn will make the algorithm slower and more unstable to fast state changes.The value ofεn 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 σs, it is possible to estimate the PC-W using the Tchebysheff’s inequality as Ωn≈6σs∝σs. Likewise the PC-G can be estimated as εn≈6σs/(M-1)∝σs. These relations show that σs is directly proportional to PC-W and PC-G. Because of the previous observation it is possible to conclude that, given a swarm of M particles, during the convergence phase a bigger value of σs would be more appropriate because it would lead to faster results. On the other hand, after the convergence phase, smaller values of σs would be more appropriate to increase the estimation precision.For these reasons in this paper a novel algorithm that manages the dynamic evolution ofσs across the iterations, namely, Particle Swarm Adaptive Scattering (PSAS), is proposed.The distance metricdn(m) 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 dn(m) is calculated between the baseline estimation α~n introduced in (5) and the resampled particle swarm αn,RES(m) as follows:(12)dn+=avgdnm+stddnm,dn-=avgdnm-stddnm,dnm=αn,RESm-α~n,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 dn- and dn+ represent, respectively, the lower and higher bounds of the error cloud dn(m) and are used as thresholds by the PSAS algorithm to trigger the increment or the decrement of σs,n as follows:(13)σs,n=σs,n-11±δσs,where δσs>0 is the proportional feedback parameter of the PSAS. The value of δσs, depending on the two feedback parameters dn+ and dn- and the threshold value γth, is set positive to iterative increment σs,n and negative to decrement σs,n and is set to zero to keep it constant.The value ofσs,n is incremented with respect to its previous value when the condition of nonconvergence is detected. This condition holds true when the error cloud dn(m) is reasonably out of [-γth,γth] convergence interval, that is, when dn->γth or dn+<-γth. Then, the PC-W is incremented to encourage convergence. If the previous condition is not satisfied, the value of σs,n is decremented w.r.t. its previous value when the convergence state is detected. This condition holds true when the error cloud dn(m) is on the edges of [-γth,γth] convergence interval, that is, when dn+>γth or dn-<-γth. Then, the PC-G is reduced to increment the precision. For the remaining values of dn+ and dn- the condition of tight convergence is reached. In this state the value of σs,n 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 δσs 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 δσs and γth 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[xn,yn,zn] at time n of the UAV evolves as(14)xn=xn-1+Tvncos⁡ϕn′sin⁡θn′,yn=yn-1+Tvnsin⁡ϕn′sin⁡θn′,zn=zn-1+Tvncos⁡θn′.The angle ϕn′, depicted in Figure 1, is the local azimuthal angle at time n of the actual trajectory coordinates [xn,yn,zn] seen from the previous path point [xn-1,yn-1,zn-1]. The local elevation angle θn′ can be defined analogously. vn is the absolute speed of the UAV. At each step, the trajectory increments generated by θn′, ϕn′, and vn are related to the “seed” components θseed,n, ϕseed,n, and vseed,n and to the random components δθn, δϕn, and δvn as follows:(15)θn′=θseed,n+δθn,ϕn′=ϕseed,n+δϕn,vn=vseed,n+δvn.The random components are generated following, respectively, the normal distributions δθ~N(mδθ,σδθ), δϕ~N(mδϕ,σδϕ), δv~N(mδv,σδv). The parameters of the random component distributions affect the UAV path behaviour. For example, for mδv>0, the UAV will have an accelerating trend and for mδϕ>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(16)θseed,n=arccos⁡ΔznNang,ΔrnNang,ϕseed,n=arctan⁡ΔynNang,ΔxnNang,vseed,n=∑k=1Nspvn-kNsp,where arctan⁡(·,·) is the arctangent function defined within [0,2π) and where(17)ΔxnNang=xn-1-xn-Nang,ΔynNang=yn-1-yn-Nang,ΔznNang=zn-1-zn-Nang,ΔrnNang=ΔxnNang2+ΔynNang2+ΔznNang2.In Figure 1, the relations among the CARATE components for the variable ϕn are represented. As shown in (16), the smaller the window width of the seed parameters Nang and Nsp, 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).Figure 1
Diagram of CARATE components for variableϕ (Nang=5).Figure 2
Example of UAV 3D trajectory generated with CARATE.The CARATE algorithm, as visible in the sample path in Figure2, generates a 3D trajectory inside a limited flight area, formed by the horizontal bound region Sxy and the vertical bound region Sz. Sxy is defined as {|x|<Lxy∩|y|<Lxy} while Sz is defined as {z>Lzmin∩z<Lzmax}. 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 xn or yn exceed Sxy the value of ϕseed,n introduced in (16) is forced for the subsequent Nbend steps to follow a set of equispaced values between ϕseed,n and ϕseed,END. The value of ϕseed,END is appropriately defined depending on the Cartesian quadrant where the trajectory crosses the surface of Sxy. Analogously, if zn exceeds Sz the value of θseed,n introduced in (16) is forced for the next Nbend steps to a set of equally spaced values between θseed,n and θseed,END. The value of θseed,END is also defined depending on the Cartesian quadrant where the surface of Sz has been crossed.


## 7. 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 [θn,ϕn] are carried out from the information given by the measurements un(x), un(y), and un(z).Now, we definebn=[bn(a)]n∈[x,y,z] as the vector of ideal measured values, namely, un(a), in (4) but without the detrimental effects due to the impairments introduced in (2) such as PN, DS, and C-AWGN:(18)bna=A2ejψ^na,where, as introduced in Section 2, ψ^n(a) is the D-PoA for ath 3A-ULA antenna branch and |A| is the constant amplitude of the complex downconverted constellation. The vector of real measurements un=[un(x),un(y),un(z)], defined in (4), can be written as(19)una=2Na∑i=1Na/2A2ejψ^na+Wna,i=A^ejψ^na+W^na,where the intermediate signal |A|2ejψ^n(a)+Wn(a,i) is averaged through antenna pairs and A^=(2/Na)∑i=1Na/2‍|A|2=|A|2. The noise components in (19) are(20)Wna,i=Aejψna,2i-1wna,2i∗+A∗e-jψna,2iwna,2i-1+wna,2i-1wna,2i∗,W^na=2Na∑i=1Na/2Wna,i,where the notation (·)∗ denotes the complex conjugate operator. wn(a,i) is the C-AWGN with zero mean and variance σw2. ψn(a,i) is the PoA as defined in (3). The expectation of the differential noise Wn(a,i) among the antennas of the same branch is mW=0, while its variance σW2=2|A|2σw2+σw4. Analogously, the average differential noise W^n(a) is described by its mean mW^=0 and its variance is σW^2=(2/Na)σW2.To calculate the CRLB it is necessary to evaluate the log-likelihoodL(ϕ,θ)=ln⁡P[un∣bn] 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 un 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(21)Lϕ,θ=-3ln⁡πσW^2-1σW^2b-u2,where the time index n is omitted until the end of the CRLB derivation for notational simplicity. The CRLB of the variable αr, defined as rth element of the unknowns vector α=[ϕ,θ], is(22)CRLBαrϕ,θ=Fϕ,θ-1r,r,where F(ϕ,θ) is the Fischer information matrix:(23)Fϕ,θ=-E∂2∂αh∂αkLϕ,θh,k∈1,2.For the log-likelihood function described in (21), we obtain(24)Fϕ,θ=∂β∂ϕ2∂β∂ϕ∂βT∂θ∂β∂θ∂βT∂ϕ∂β∂θ2=22πA2DnormσW^2sin⁡θ2001,where β=[cos⁡ϕsin⁡θ,sin⁡ϕsin⁡θ,cos⁡θ] is the unit direction vector and Dnorm=D/λ0 is the normalized interantenna distance (Figure 4). The notation (·)T denotes the transpose operator. Finally, the CRLB expressions for ϕ and θ are(25)CRLBϕϕ,θ=12σW^2πA2Dnorm21sin⁡θ2,CRLBθϕ,θ=12σW^2πA2Dnorm2.The RMSE of every estimator of the unknown αr from the measurements u is lower bounded by CRLBαr. 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 un(a) introduced in (4) at time nT, 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.


## 8. Performance Analysis

The performance of the PF algorithm is now assessed. A comparison with the BM estimation introduced in Section3 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 N time samples long. In order to have a common performance evaluation metric, in the following we define the average of the Root Mean Square Error (RMSE) calculated through every trajectory sample, namely, the Average Trajectory RMSE (AT-RMSE).n t r a jth generated trajectory is described with the set {θn(ntraj),ϕn(ntraj)}, its BM estimation is {θ~n(ntraj),ϕ~n(ntraj)}, and its PF estimation {θ^n(ntraj),ϕ^n(ntraj)}. Their AT-RMSE are, respectively, RMSE¯α,BM=(1/Ntraj)∑ntraj=1Ntraj‍RMSEα,BM(ntraj) and RMSE¯α,PF=(1/Ntraj)∑ntraj=1Ntraj‍RMSEα,PF(ntraj) for α∈[θ,ϕ] where(26)RMSEα,BMntraj=∑n=1Nαnntraj-α~nntraj2N,RMSEα,PFntraj=∑n=1Nαnntraj-α^nntraj2Nare the RMSEs calculated for each trajectory ntraj. Ntraj 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), asCRLB¯α=(1/Ntraj)∑ntraj=1Ntraj‍CRLBα(ntraj) for α∈[θ,ϕ] and where(27)CRLBαntraj=1N∑n=1NCRLBαϕnntraj,θnntraj.The default parameters for the performance evaluation are (a) anglesσδθ=7°, σδϕ=7°, and Nang=20; (b) speeds vmin=0.5 km/h, vmax=10 km/h, and Nsp=2. The array parameters are Na=6 and Dnorm=1/3. The PF algorithm basic settings are M=20 and Nmm=20. PSAS parameters δσs=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=[|A|2/σw2] dB =10 dB. The number of simulated trajectories Ntraj and their length N have been properly dimensioned to guarantee the statistical confidence of the results.In Figure2 an example of 3D-trajectory generated for performance evaluation is shown. The 3A-ULA is positioned in the center of the axis.


### 8.1. Performance with Static and Dynamicσs

We now consider the performance as a function of the signal-to-noise ratioSNR=|A|2/σw2. 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 σs, 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 σs. 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.Figure 3
AT-RMSE for differentSNRdB conditions with dynamic and static σs.Figure 4
AT-RMSE for different normalized 3A-ULA interantenna distanceDnorm=D/λ0.


### 8.2. Performance for Different Interantenna DistancesDnorm

Now, the behaviour of the AT-RMSE for the BM and PF algorithms is analysed varying the normalized interantenna distanceDnorm. To unambiguously extract the phase from un(a) the interelement distance must be less than half of the impinging signal wavelength [15]: Dnorm≤0.5. Furthermore in (25) it is visible that for higher values of Dnorm the CRLBs increase. This leads to Dnorm=0.5 as the optimal value for the estimation. However, in lower SNR scenarios and for Dnorm 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 un(x) 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 Dnorm=Df0/c0 than the BM. Thus, the same 3A-ULA is usable for a wider set of carrier frequency f0 with respect to the BM that exhibits its optimal performance for a narrow interval of Dnorm 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.


## 8.1. Performance with Static and Dynamicσs

We now consider the performance as a function of the signal-to-noise ratioSNR=|A|2/σw2. 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 σs, 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 σs. 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.Figure 3
AT-RMSE for differentSNRdB conditions with dynamic and static σs.Figure 4
AT-RMSE for different normalized 3A-ULA interantenna distanceDnorm=D/λ0.


## 8.2. Performance for Different Interantenna DistancesDnorm

Now, the behaviour of the AT-RMSE for the BM and PF algorithms is analysed varying the normalized interantenna distanceDnorm. To unambiguously extract the phase from un(a) the interelement distance must be less than half of the impinging signal wavelength [15]: Dnorm≤0.5. Furthermore in (25) it is visible that for higher values of Dnorm the CRLBs increase. This leads to Dnorm=0.5 as the optimal value for the estimation. However, in lower SNR scenarios and for Dnorm 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 un(x) 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 Dnorm=Df0/c0 than the BM. Thus, the same 3A-ULA is usable for a wider set of carrier frequency f0 with respect to the BM that exhibits its optimal performance for a narrow interval of Dnorm 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.


## 9. 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 the integration of ranging algorithms like Received Signal Strength Estimation (RSSE) and ToA to PF to provide full 3D localization.



---

*Source: 7630950-2016-09-21.xml*