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

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 . 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 introduced in (16) is forced for the subsequent steps to follow a set of equispaced values between and . The value of is appropriately defined depending on the Cartesian quadrant where the trajectory crosses the surface of . Analogously, if exceeds the value of introduced in (16) is forced for the next steps to a set of equally spaced values between and . The value of is also defined depending on the Cartesian quadrant where the surface of 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 are carried out from the information given by the measurements , , and .

Now, we define 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 , defined in (4), can be written aswhere the intermediate signal is averaged through antenna pairs and . The noise components in (19) arewhere the notation denotes the complex conjugate operator. is the C-AWGN with zero mean and variance . is the PoA as defined in (3). The expectation of the differential noise among the antennas of the same branch is , while its variance . Analogously, the average differential noise is described by its mean and its variance is .

To calculate the CRLB it is necessary to evaluate the log-likelihood of the acquired data given the parameters to be estimated. represents the natural logarithm. In this context, for the calculation of CRLB, the measurements vector 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 thatwhere 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 , iswhere is the Fischer information matrix:For the log-likelihood function described in (21), we obtainwhere is the unit direction vector and is the normalized interantenna distance (Figure 4). The notation denotes the transpose operator. Finally, the CRLB expressions for and areThe RMSE of every estimator of the unknown from the measurements is lower bounded by . 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 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.

8. 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 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).

th generated trajectory is described with the set , its BM estimation is , and its PF estimation . Their AT-RMSE are, respectively, and for whereare the RMSEs calculated for each trajectory . 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), as for and where

The default parameters for the performance evaluation are (a) angles , , and ; (b) speeds  km/h,  km/h, and . The array parameters are and . The PF algorithm basic settings are and . PSAS parameters and were numerically optimized for the CARATE algorithm. In the signal model introduced in (2), the PN is considered white with standard deviation , the PO is considered uniform in for each array element, and the signal-to-noise ratio is  dB =10 dB. The number of simulated trajectories 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.

8.1. Performance with Static and Dynamic

We now consider the performance as a function of the signal-to-noise ratio . 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 for and 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 for and 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.

8.2. Performance for Different Interantenna Distances

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

Competing Interests

The authors declare that they have no competing interests.