Multi Direction-of-Arrival Tracking Using Rigid and Flexible Antenna Arrays

This paper is concerned with the problem of simultaneously tracking the direction-of-arrival (DOA) of far-field multiple moving sources/users in wireless communications using the vector-signal received by an antenna array of $N$ elements. The antenna array can be rigid (fixed array locations) or flexible (time-varying array locations), and it is used in conjunction with a “manifold extender”, a spatiotemporal state-space model and a Kalman-type tracking approach for non-stationary wireless channels. In particular, two tracking approaches are proposed. The first is based on an arrayed Extended Kalman Filter (arrayed-EKF) algorithm and the second on an arrayed Unscented Kalman Filter (arrayed-UKF) algorithm. Furthermore, if the array is rigid the spatiotemporal state-space model incorporates the DOAs and the angular velocities of the sources, while if it is flexible it also includes the array locations in the set of state-variables. The performance of the two approaches using both rigid and flexible arrays is evaluated using computer simulation studies and compared with a subspace tracking algorithm and a particle filter method under the same conditions. The results show that the arrayed-UKF and the arrayed-EKF show superior tracking performance, especially for low SNRs.

ing sources/users using an array has been an important research area with a wide range of applications in sonar, radar, air traffic control, wireless communications, remote sensing, etc.. DOA tracking techniques can be classified into probabilistic and parametric. For instance, in [1], a probabilistic DOA tracking approach is presented which is based on a probability hypothesis density filter with a likelihood function expressed as a complex Wishart random matrix. In [2], the tracking of DOAs is also probabilistic and is based on the sparse approximation technique LASSO. One of the main tracking families of parametric techniques is the "subspace tracking" [3]- [6]. In [3], a cross-correlation based 2D DOA tracking algorithm is proposed using an automatic pair-matching batch method which, however, is restricted to only L-shape array geometries. Some "subspace tracking" techniques are based on various decomposition forms such as Singular Value Decomposition (SVD) [4], URV decomposition [5] and cross RV decomposition (CRV) [6].
However, many DOA tracking techniques assume that the sources are stationary over a small time frame (observation interval) and for each time frame a DOA subspace estimation algorithm like multiple signal classification (MUSIC) [7] can be applied. In these type of techniques the tracking is based on repetitive DOA estimation (or, in general, repetitive localization) where the set of consecutive estimates provide the tracking trajectory (e.g. [8]). However, in non-stationary environments, repetitive high-resolution estimation algorithms for DOA source trajectory tracking exhibit serious performance degradation. In addition, these techniques suffer from the data association problem and several algorithms have been proposed to avoid this problem. In [9], the authors minimize the norm of an error matrix based on the covariance matrix of the received array output. A repetitive DOA source tracking algorithm is proposed in [10] which uses the most recent received data to update the existing DOA estimates using the MUSIC algorithm. Both [9] and [10] avoid data association by preserving the order of the estimated DOAs over certain iterations, which however suffer from spread array spatial spectrum This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ effects caused by the source motion. Alternatively, statespace model based approaches have been proposed which are combined with various tracking algorithms (e.g. [11]- [14]). In [11], multiple target states (MTS) are defined to describe the source motion, followed by a ML algorithm for updating the MTS and tracking the DOAs. In [12], a bank of linear combiner matrices are formed as state variables, and then they are updated adaptively by an H-infinity filter to track the noise subspaces and consequently the DOAs. Furthermore, some Bayesian state-space models have been used in [13], [14] which incorporate the source motion and the likelihood of received measurements based on the source's state. Then, a particle filter is used to track the source state to obtain the DOA of the sources. The particle filter has also been used in [15] and [16] for target tracking in radar system. However, in radar literature, "active" radar tracking approaches have been employed where both system architectures and assumptions are different than those used in communications.
Another family of approaches for repetitive DOA estimation are based on Kalman Filter (KF). Since the standard Kalman filter requires the measurement model to be linear, all the approaches that use standard Kalman filter need to pre-estimate the DOA and view the DOA estimates as "measurement". Then the pre-estimated DOA is refined by the Kalman filter. For example, the Kalman filter has been used for multiple DOA tracking which are pre-estimated by a least squares (LS) estimator [17], or a Maximum Likelihood (ML) estimator [18], [19]. In [20], the proposed algorithm improves the algorithm in [10] by employing a source movement model and a Kalman filter. The Kalman filter has also been used in [21] for tracking signal subspace towards the objective of tracking single target. However, it is important to point out that the above KF approaches require the overall observation interval to be divided into small intervals over which the DOAs can be assumed to be stationary. Consequently, they will suffer from serious performance degradation when this assumption is not valid.
In addition to KF, the Extended Kalman filter (EKF) and the Unscented Kalman Filter (UKF) are suitable for the case that the measurement model is non-linear. For instance, the EKF has been employed in [22] for trajectory tracking of moving sources using a large aperture rigid array and in [23] for DOA tracking of moving sources using a small aperture rigid array. In [24], the EKF is combined with a particle filter forming an EKPF algorithm. However, to the best of our knowledge, there are not many papers for DOA tracking using EKF and the majority of these papers are for single DOA tracking. In addition, the UKF has been rarely used and in [25], [26] it has been used for single DOA tracking. In general, the EKF and UKF algorithms have equal computational complexity but they are conceptually different. The EKF linearizes nonlinear transformations by using Taylor series expansions and then uses these linear transformations in standard Kalman filter. Whereas the UKF involves the unscented transformation which essentially selects a set of points (sigma points) via a deterministic sampling approach, and then propagates these points to the true nonlinear function which are then exploited to form the mean and covariance of the estimation [27], [28]. In this paper, both the EKF and UKF are employed to track multiple DOAs in non-stationary environment and our paper is the first paper to employ UKF for multi-DOA tracking.
All the above methods, like the majority of array processing techniques assume by default that the array geometry is rigid. On the other hand, using flexible array geometries is an interest problem in array signal processing for airborne, vehicular, underwater and other applications [29]. The "flexible array" is defined as an antenna array with time varying geometry, i.e. each of the array elements move independently and this paper is an extension of [30] where a flexible array formed by a swarm of unmanned aerial vehicles (UAVs) is presented. If the array geometry is flexible, i.e. the array geometry changes as a function of time, then the majority (if not all) of the array processing algorithms and theory cannot be directly used.
Several approaches have been proposed for solving source tracking problem with time varying arrays but for static sources. For example, the ML estimator has been employed in [31]- [33]. However, in these cases, a multi-dimensional search is required which is computationally prohibitive. In [34], two eigenstructure-based algorithms based on the concept of array interpolation and focusing matrices are proposed with faster approximations to the ML estimators. However, these algorithms need relatively large Signal to Noise Ratio (SNR) levels to maintain satisfactory performance. In [35], the authors use noncoherent time-varying arrays that are a collection of coherent subarrays with stationary covariance matrices. Then, the functions of the covariance matrices of the subarrays are derived to find source locations. As it was stated before, all these approaches ([31]- [35]) are designed for static sources. In [36], a dynamic radar network is proposed which is related with a swarm of UAVs working independently for tracking the location of a single target. It uses a Markovian state space model and employs a two-step EKF algorithm. However, this approach is not related with flexible arrays as the UAVs work independently and this scenario belongs to radar applications. This paper is concerned with tracking the DOA of multiple sources in non-stationary wireless communications using both rigid and flexible arrays. In addition to [30] in the current literature, to the best of our knowledge, there are only two papers [37], [38] that deal with the tracking using timevarying arrays. Reference [37] is concerned with tracking stationary sources using an array of hydrophones. This is a "towed-array" -towed behind a submarine or a surface ship on a cable -where the hydrophones are placed at specific/constant distances along the cable. This is also a flexible array because when the ship turns this line becomes curvy and there are small but very restricted changes in the overall shape. However, [37] deals with DOA estimation of stationary acoustic sources using Maximum Likelihood (ML) followed by a second estimation of the shape of the array using Kalman filtering. Reference [38] is a "probabilistic" approach for non-stationary channels which completely ignores the parameterisation in terms of the array manifold vector and array geometry. It recursively estimates a conditional probability density function (Bayesian filtering) by following the EKF iterative steps, although this is not an EKF approach, while its performance is compared with that of a particle-filter method.
In our proposed approaches, all the sources and the array locations are tracked in a unified way which is suitable for tracking even of fast moving sources using antenna array systems. In particular, two novel approaches are proposed based on an arrayed-EKF and an arrayed-UKF using • a rigid array and a flexible array, • a spatiotemporal state-space model and • a "manifold extender" for simultaneously tracking multiple DOAs snapshot-bysnapshot. In our proposed approaches, both "rigid" and "flexible" arrays as well as arrayed-EKF and arrayed-UKF are presented in a unified way. Furthermore, the concept of "manifold extender" (see [39]) is employed which increases the "degrees of freedom" of the system by integrating the spatial and the temporal domains. This is used in the spatiotemporal state-space model. Note that the "extended manifold" vectors are longer vectors than the spatial manifold vector and the locus of these vectors are mathematical objects (manifolds) that are embedded in a larger observation space than the one provided by the number of antenna array elements N (see [39]). If the array is flexible, apart from tracking multiple DOAs, the array locations are also simultaneously tracked as they change arbitrarily with time. Therefore, we developed our array processing algorithms (arrayed-EKF/UKF), based on fixed or time-varying array geometry, where all the antenna array elements (in a constant or a time varying geometry) work together as one unit for solving the problem of trajectory tracking of moving sources. The integration of all the above forms our proposed "arrayed-EKF" and "arrayed-UKF" algorithms.
Note that the proposed algorithms in this paper have many applications, including UAV communications. For instance, a rigid antenna array may be deployed on a single UAV platform 1 for tracking multiple users. Furthermore, a number of UAVs (a swarm of UAVs), each equipped with a single antenna having its own propulsion system, can be used as a flexible array for multi-user tracking (e.g. [30]). In this case, each UAV is equipped with a GPS-clock so that the array system will have a common clock to keep the system coherent.
The remainder of this paper is organized as follows. In Section II, the array system and the received vector-signal model are presented. In Section III, the mobility model of the multiple sources is described. Then, the extended mobility model is presented. In the case of flexible arrays, the mobility model of the array elements is also discussed. In Section IV, the proposed approaches are introduced based on arrayed-EKF and arrayed-UKF for both rigid and flexible arrays. In Section V, the performance of the proposed approaches is evaluated using computer simulation studies. Finally, this paper is concluded in Section VI. II. SYSTEM MODEL Consider an array system for tracking multiple users with a receiver array of N antennas. Figure 1 shows the baseband representation of the array system model consisting of M far-field transmitters/users, a wireless channel and an array receiver. 2 With reference to Figure 1, at Point-A, the transmitted data sequence of symbols of the i-th user {a i [n]} (with a symbol duration of T cs ) where each symbol is weighted by a N c × 1 weight-code vector given by Then, at Point-B, the weighted data symbol sequence } is driven to a Digital-to-Analog Converter (DAC) to produce a baseband transmitted signal m i (t) at Point-C. At the receiver, an antenna array of N antennas is employed with locations where the vector r m ∈ R 3×1 denotes the Cartesian coordinates of the m-th antenna and the N × 1 vectors r x , r y , r z are the Cartesian coordinates of all antennas on the x-axis, y-axis and z-axis, respectively. In this paper, we also consider the scenario where the array locations (and thus the array geometry) change due to any unknown forces. In this case, the array is flexible with time-varying geometry, re-modelled as functions of time as follows Thus, at Point-D in Figure 1, the received baseband vector signal x (t) ∈ C N ×1 can be modelled as follows where, for the i-th user, β i (t) denotes the path fading coefficient, F i represents the Doppler frequency and the vector n(t) ∈ C N ×1 denotes the additive white complex Gaussian noise of zero mean and covariance matrix given by with σ 2 n denoting the noise power. In Equ 5, S i (t) ∈ C N ×1 denotes the time-varying manifold vector which is given as follows for flexible array  where F c is the carrier frequency, c denotes the speed of light and u i (t) is given by with θ i (t) and φ i (t) representing the azimuth and the elevation angles of the i-th source. In general, the vector u = u (θ, φ) denotes the (3 × 1) unit-norm vector pointing towards the direction (θ, φ), as illustrated in Figure 2. In this paper, with no loss of generality, the elevation angle is assumed to be zero (i.e. φ i (t) = 0). Equ 5 can also be rewritten in a more compact form as follows where S(t) ∈ C N ×M is the matrix with columns the array manifold vectors, i.e.
and m(t) ∈ C M×1 is expressed to include, in addition to the M baseband message signals, the Doppler frequencies and the path coefficients as follows with its covariance matrix R mm defined as With reference to Figure 1, the (N × 1) vector signal x (t) is firstly discretised (see Point-E). Then, at Point-F the (N N ext × 1) vector signal x [n] can be expressed as follows where the vector n [n] ∈ C N Next×1 denotes the discretised "extended" noise (i.e. the noise at Point-F in Figure 1) and the matrix H [n] ∈ C N Next×M is given by with its i-th column h i [n] ∈ C N Next×1 denoting the timevarying extended manifold vector of the i-th user given as follows with N ext = N c . Furthermore, at the transmitter, if the user being tracked does not include the weight w i , then w i = 1 Next may be used. In this case, the extended manifold vector of Equ 15 is simplified to Therefore, the extended manifold vector here has two forms/cases as described in Equ 15 and Equ 16 although other forms from [23], [39] may be included. In addition, if the manifold extender is not employed in this model, then N ext = 1 and Thus, the vector h i [n] ∈ C N Next×1 can be expressed as follows Note that the array locations that constitute the extended manifold vector h i [n] may be a function of time (changing from symbol to symbol), depending on whether the array is rigid or flexible.

A. Source/User Mobility Model
Since the array operates in the presence of M co-channel users, Figure 3 illustrates the mobility model of the i-th user, for i = 1, 2, . . . , M relative to the array reference point. As shown in Figure 3, the i-th user moves in an arbitrary direction with velocity-vector v i ∈ R 3×1 given as follows where the velocity vector v i is decomposed into two other vectors: • the radial component v ρi u vρ i , and • the orthoradial (angular) component v θi u v θ i , with v ρi denoting the radial velocity of the i-th source and v θi representing its angular velocity. The vectors u vρ i and u v θ i are unity vectors that are mutually orthogonal. The radial velocity v ρi causes the Doppler effects which is modelled as exp(j2πF i t) in Equs 5 and 11. Note that F i is given by As shown in Equ 13, the Doppler coefficient has been incorporated into the combined vector-signal m[n] which will then be estimated and utilised by the proposed tracking algorithms. Thus, the radial velocity does not affect the DOA/azimuth θ i . If the i-th user moves with an angular velocity v θi , which can be constant or variable, its azimuth angle may be described as where is the time elapsed between t and t 0 . If T is assumed to be equal to the sampling period, i.e.
then Equ 21 is discretised and becomes Thus, let us define the discrete-time state vector for the i-th Then, the discrete time kinematic model for the i-th user for t = nT is given as where G i ∈ R 2×2 is the transition matrix given by In Equ 27, b i [n] ∈ R 2×1 represents perturbations about the azimuth angle and the angular velocity, and can be modelled as "noise" with zero mean and covariance matrix R bi bi ∈ R 2×2 given by where σ 2 θi denotes the continuous time model intensity for the azimuth-velocity.

B. Flexible Array Mobility Model
With reference to Figure 3, it is shown that the array at the receiver can be either rigid or flexible. The movement of the flexible array may include • a known motion of the whole array which is represented by the motion of the array's reference point and does not affect the array geometry, • a known motion of its individual elements relative to the reference point, and • small unknown motion or perturbations casued by any unknown forces 3 with the last two motions changing the array geometry. Thus, in a flexible array the tracking of the array geometry (i.e. the tracking of the locations of the array elements) is essential.
With a sampling period T , the discrete time mobility model of the array locations 4 is expressed as follows where the vector b xy [n] is given as follows which includes the instantaneous array locations on the x-axis and y-axis, respectively. The vector b xy [n] ∈ R 2N ×1 in Equ 30 denotes the perturbations associated with their respective locations, with its intensity σ 2 xy , and its covariance matrix is given as follows The matrix G xy ∈ R 2N ×2N is a block diagonal matrix containing known transition matrices of all the array elements shown as where F j ∈ R 2×2 represents a known transition matrix of j-th array element. For instance, if ω j denotes the angular velocity of the j-th array element about the reference point, then

C. Overall Mobility Model
The overall mobility model, which is the discrete time M -user kinematic model, is constructed based on Equs. 27 and 30, as follows where b[n] is the overall discrete-time state vector which is constructed as T ∈ R (2N +2M )×1 flexible array (36) and its perturbation vector b[n] is given as follows 3 Note that, even in fixed array geometries, constant uncertainties in the array locations may decrease the performance of the direction-finding system [31]. 4 With no loss of generality, the antenna array elements are assumed to be located on the (x, y) plane, i.e. r z [n] = 0 N .
where N dim = 2M rigid array 2N + 2M flexible array (40) and the perturbation matrix R b b is given as Table I provides the dimensionality of the various vector/ matrix parameters used in the mobility models. The mobility models employed in this section can be seen as some representative examples but other mobility models may be used in this proposed framework.

IV. TRACKING ALGORITHMS BASED ON SPATIOTEMPORAL STATE-SPACE MODEL
In this section, based on the antenna array models presented in Sections II and III, two tracking algorithms are proposed which belong to the Kalman family of techniques for non-linear "measurement" models. We will call them "arrayed-EKF" and "arrayed-UKF" as they are based on the integration of • the array signal model of Section II for both "rigid" and "flexible" antenna arrays, • the mobility models of Section III, and • the EKF/UKF iterative theoretical tools. Based on the snapshot at Point-F in Figure 1, modelled by Equ 13, and the overall mobility model given by Equ 36, the spatiotemporal state-space model is constructed as follows describing the dynamics of all users' motion (Equ 42) and the vector signal received by the rigid or flexible arrays (Equ 43). The two proposed algorithms are summarised in Tables II and III. Note that the notation has been simplified, by replacing the symbol index [n] with a subscript n and, thus, the following notation is employed in Tables II and III It is important to point out that, in the presentation of the two algorithms in these two tables, the selection of a common Kalman structure is deliberately maintained for better clarification of each step. However, in the arrayed-UKF, which is based on the square-root UKF, the state-vector b [n] of Equ 36 and the non-linear "measurement" vector x [n] of Equ 43 should be further processed to form two matrices B n and X. In particular, the matrix B n ∈ R N dim ×(2N dim +1) can be formed as follows where • B j is the j-th column of the matrix B n which is known as the j-th "sigma point-vector". • η is a scaling factor • T n is an N dim × 2N dim matrix which is the Cholesky factorisation 5 of the covariance matrix P n of the state vector of Equ 36. That is Then, by applying the nonlinear function to these "sigma point vectors" the measurement matrix X ∈ C N Next×(2N dim +1) is formed, i.e.
Thus, using the unscented transformation instead of the Jacobian matrix of the nonlinear measurement, and propagating the 5 The Cholesky factorisation of Pn is unique as Pn is a positive definite matrix.
Cholesky factor T n instead of the covariance of the estimate error P n , the proposed arrayed-UKF algorithm is summarised in Table III. Note that in Table III, the matrices W (m) and W (c) are diagonal matrices whose diagonal values are the weights to compute the mean and the covariance of the measurement, respectively. These matrices have the following definitions: where the scaling parameters μ and β are as follows with the constant α controling the spread of the "sigma pointvectors" around b n , and ρ compensating for the distribution of b n . Furthermore, the parameter η in Equ 46 is related to μ as follows With reference to Table II and Table III, the initial state b 0 can be provided by a priori guess or pre-estimated by any kind of the DOA estimation approaches, such as the ML, LS, or subspace approaches. The initial covariance P 0 can be set to γI 2M where γ indicates the confidence level in the accuracy of the initial estimates.

V. COMPUTER SIMULATION STUDIES
The performance of the proposed algorithms are evaluated in this section using computer simulation studies. The data symbol sequence transmitted by each user is assumed to be a random complex sequence of zero mean and unity variance. The weight-codes are gold-codes of ±1s of length N ext = N c generated by modulo-two addition of two m-sequences described by the polynomial D 3 +D 2 +1 and D 3 +D+1. This implies that N ext = N c = 7. The user tracking is assumed to be carried out over a time interval of 5000 spatiotemporal snapshots x [n] ∈ C N Next×1 collected at Point-F in Figure 1. The continuous time model intensity for the azimuth-velocity (see Equ 29) is set to σ 2 θi = 1.1 × 10 −8 (deg /T ) 2 , ∀i and the intensity of the perturbations for the array locations is set to σ 2 xy = 6.4 × 10 −7 (λ/T ) 2 . The parameter γ of the initialization stage of the proposed algorithms is set to 10 −6 . The parameters used in Equs 51 and 52 are set to α = 10 −4 , and ρ = 2.

A. Rigid Array Geometry
For the rigid antenna array, the geometry is assumed (without any loss of generality) to be a grid planar array of 9 elements and its locations are shown in Figure 4. Furthermore, it is assumed that the array operates in the presence of 4 far-field moving sources/users and their initial angular velocities are assumed to be 0 deg /T which then change according to the velocity trajectory shown in Figure 5. Figure 6 shows an example of DOA trajectory tracking of four far field moving sources using the proposed arrayed-EKF and arrayed-UKF approaches for SNR = 10 dB. In Figure 6, the tracking results of the "Source 2" between the 3700-th snapshot and the 3800-th snapshot (framed area) are zoomed to show the tracking performance using both approaches. It is clear that both approaches track the DOAs with high accuracy, especially the arrayed-UKF approach. Then, the performance of the proposed arrayed-EKF and arrayed-UKF approaches for all three cases of Equ 18 is examined using Monte Carlo simulation studies. The results are shown in Figure 7 where the Root Mean Square Error (RMSE) of the estimated azimuth angles is plotted as a function of the SNR for 500 Monte Carlo simulations, where in each simulation the error is averaged over the whole trajectory corresponding to 3000 snapshots. It is evident in Figure 7  • the subspace tracking algorithm presented in [12] which employs an H-infinity filter for estimating the noise sub- space which is then used to track snapshot-by-snapshot the DOAs of multipaths, • the particle filter approach of [40], but integrated with the received array vector-signal model for both rigid and flexible antenna arrays and the overall mobility model with the state vector of Equ 36. This particle filter will be referred here as "arrayed-PF". The comparative results are shown in Figure 8 and it is evident that the proposed algorithms outperform both the subspace tracking and the "arrayed-PF"algorithms.
Next the RMSE performance of the arrayed-UKF algorithm for the extended manifold vector h i [n] = S i [n] ⊗ w i is examined. Figure 9 shows the RMSE of the estimated DOA angles as a function of the number of spatiotemporal snapshots using the proposed arrayed-UKF algorithm for different SNRs. Note that, in Figure 9a, the initial DOAs are assumed known, and in Figure 9b the initial DOAs are obtained from a random Gausian distribution with unity variance. The results show that the estimation error for different signal environments is small and remains constant over time, which illustrates the stability of the proposed arrayed-UKF approach. Overall, based on    RMSE of the estimated source azimuth angles averaged over 3000 spatiotemporal snapshot evaluations using the arrayed-EKF and the arrayed-UKF approaches for the three cases of Equ 18 (500 iterations). Fig. 8. Comparison of the proposed arrayed-EKF and arrayed-UKF algorithms with other algorithms (the subspace tracking [12] and an arrayed-PF algorithm based on [40] using 100 particles). the results of both Figures 7 and 9, it can be concluded that the arrayed-UKF algorithm offers significant tracking accuracy over the SNR range, which indicates its robustness to noise.

B. Flexible Array Geometry
For the flexible antenna array, the initial geometry is assumed to be also given by Figure 4 but then the geometry will change randomly and a representative trajectory of the 9 antennas relative to the array reference point is shown in Figure 10. The fifth antenna array element is assumed to be the array reference 6 point. For this trajectory, the angular velocity for each flexible array element is assumed to be ω j = 0.01(deg /T ), ∀j. The angular velocities of the sources are shown in Figure 11. Figure 12 shows another example of DOA trajectory tracking of four moving sources using the proposed arrayed-EKF and arrayed-UKF approaches with a flexible array geometry for SNR = 10 dB. The tracking results of "Source 1" between the 3700-th snapshot and the 3800-th snapshot are zoomed in. It can be observed that the DOA tracking works well with a flexible array. The estimated array locations for n = 3000 and 5000 by the two algorithms are illustrated in Figures 13 and 14, respectively. Furthermore, the trajectories of the flexible array locations and the trajectories of the estimated array locations using the arrayed-UKF algorithm are shown in Figure 15. It is clear that the array geometry changes dramatically with time and its instantaneous array locations are successfully estimated by the proposed approaches.
Finally, for the flexible array case, 500 Monte-Carlo simulations have been carried out under the same simulation environment as used in Figure 12. The two proposed algorithms are also compared with the "arrayed-PF" approach (using 150 particles) for the flexible array case under the same simulation environment. Note that the "subspace tracking" algorithm [12] examined in Figure 8 is not able to work 6 The array reference point also moves but this motion is not shown (only the relative motion with respect to this reference point is shown).
in the case of flexible arrays. The results are shown in Figures 16 and 17 where the RMSE of the estimated array locations ( Figure 16) and azimuth angles ( Figure 17) are plotted as a function of the SNR where in each simulation the error over the DOA trajectory of 3000 snapshots is averaged for each moving source. These figures indicate that the proposed arrayed-UKF algorithm with the extended manifold vector h i [n] = S i [n] ⊗ w i has superior tracking accuracy in both the estimated array locations and the estimated azimuth angles.
Finally, it can be seen from Figures 6 and 12 that the proposed algorithms do not suffer from the data association problem when the source trajectories cross each other. This is because of the use of the spatiotemporal state space model which provides a one-to-one mapping between the sources and their manifold vectors.

VI. CONCLUSION
In this paper, a theoretical framework is presented for tracking far-field sources in an non-stationary environment using both rigid and flexible antenna array geometries. The proposed multi-source tracking framework is based on the integration of a spatiotemporal state-space modelling, the extended manifold concept and EKF/UKF theoretical iterative tools, for both rigid and flexible array geometries. The performance of the proposed approaches was examined using computer simulation studies under various noise levels and compared with a subspace tracking and a particle filter algorithms. The results indicate that the arrayed-UKF algorithm has better tracking performance than the other algorithms for both rigid and flexible array geometries, while it shows robustness to noisy environments.
For further extension of the current work, antenna arrays may be deployed at each transmitter of Figure 1 (i.e. MIMO). In such system, the manifold vector S i [n] in Equ 18 can be replaced by S * i [n] ⊗ S i [n], which is known as "virtual array manifold", where S i [n] denotes the manifold vector of the transmitter. In this case, the tracking performance could be further improved.    Array geometry at n = 3000 showing the initial array locations (circle), the true array locations (square) and the estimated array locations using the arrayed-UKF (diamond) and the arrayed-EKF (cross).    16. RMSE of the estimated array locations averaged over 3000 spatiotemporal snapshot evaluations using the arrayed-EKF approach and the arrayed-UKF approach for the three cases of Equ 18, and the "arrayed-PF" approach, as a function of SNR (500 iterations). Fig. 17. RMSE of the estimated source azimuth angles with the flexible array averaged over 3000 spatiotemporal snapshot evaluations using the arrayed-EKF approach and the arrayed-UKF approach for the three cases of Equ 18, and the "arrayed-PF" approach, as a function of SNR (500 iterations).