A Sensor-Driven Analysis of Distributed Direction Finding Systems Based on UAV Swarms

This paper reports on the research of factors that impact the accuracy and efficiency of an unmanned aerial vehicle (UAV) based radio frequency (RF) and microwave data collection system. The swarming UAVs (agents) can be utilized to create micro-UAV swarm-based (MUSB) aperiodic antenna arrays that reduce angle ambiguity and improve convergence in sub-space direction-of-arrival (DOA) techniques. A mathematical data model is addressed in this paper to demonstrate fundamental properties of MUSB antenna arrays and study the performance of the data collection system framework. The Cramer–Rao bound (CRB) associated with two-dimensional (2D) DOAs of sources in the presence of sensor gain and phase coefficient is derived. The single-source case is studied in detail. The vector-space of emitters is exploited and the iterative-MUSIC (multiple signal classification) algorithm is created to estimate 2D DOAs of emitters. Numerical examples and practical measurements are provided to demonstrate the feasibility of the proposed MUSB data collection system framework using iterative-MUSIC algorithm and benchmark theoretical expectations.


Introduction
DOA estimation of source using array sensors plays an important role in the field of array signal processing. DOA estimation includes one-dimensional (1D) DOA estimation (azimuth) and 2D DOA estimation (azimuth and elevation). The accuracy of DOA estimation is mainly impacted by the algorithm, the geometry of sensor array structure, signal-to-noise (SNR), and snapshots, etc. Algorithms and array geometry are two essential research topics.
The DOA estimation method has been studied for decades and it is still an active research topic in recent years. Initially, DOA estimation based on sensor array structures used the Bartlett beamforming method, but it cannot offer high-resolution due to the Rayleigh limit [1,2]. Then, Burg came up with the maximum entropy (ME) method, which is a high-resolution method, but it has a low robustness and a considerable computation [3]. Later, a series of high-resolution spatial spectrum estimation methods based on decomposition of matrix eigenvectors for array signal processing came out and created a new era. All those spatial spectrum estimation methods were represented by MUSIC and estimation signal parameters via rotational invariance technique (ESPRIT) [4][5][6][7]. These two methods have greater resolution and accuracy than other classical methods. The simulations in [8] indicate that MUSIC algorithm is more accurate and stable than ESPRIT algorithm for uniform linear array (ULA). Furthermore, ESPRIT can only be used for invariant geometry, while MUSIC can be applied for arbitrary non-uniform sensor arrays and multiple-source estimation. The MUSIC algorithm attracted more and more attention since it came out, which was the milestone of spatial spectrum estimation methods. In certain conditions, MUSIC algorithm is 1D implementation of maximum likelihood (ML) method, which shares the same characteristic with ML [9,10]. Then, due to the relatively high calculation complexity of MUSIC, some search-free algorithms, such as root-MUSIC, manifold separation based on root-MUSIC, use the root-solving technique to reduce the computational complexity [11,12]. Nevertheless, it can only be used for ULA. In order to deal with arbitrary non-uniform arrays, the Fourier domain root-MUSIC (FD root-MUSIC) algorithm was addressed [13], but the FD root-MUSIC algorithm is mainly used for 1D DOA estimation.
Most of algorithms above are still in simulation and theoretical stages. Recently, practical implementation of the DOA estimation system using field programmable gate array (FPGA) and digital signal processing (DSP) are reported [14][15][16], those implementations achieve real-time application of DOA estimation based on the MUSIC algorithm. Different from other well-known subspace-based techniques like ESPRIT, MUSIC algorithm has many advantages in the real implementation owing to its simplicity and suitability for parallel processing [16].
Apart from DOA estimation methods, the geometry of the antenna array is also very important. Most of the investigations on the array structure are limited in either 1D linear array or 2D planar array configurations [17]. The 1D array can only detect 1D DOA, and 2D array structures such as uniform rectangular array (URA) and uniform circular array (UCA) can well estimate azimuth angles but cannot well estimate elevation angles due to its small antenna aperture in the elevation direction. In order to improve the elevation angle estimation accuracy, one may develop a three-dimension (3D) array structure by putting more elements in the vertical direction to make a large elevation aperture [18,19]. However, it requires additional hardware and computational cost. Furthermore, the uniform linear and planar array will cause the ambiguity problem (angle aliasing) due to the symmetry of array structures [20][21][22]. Xia et al. proposed that cubic arrays still have ambiguity problem and the spherical array can significantly reduce this problem [22]. Recently, 3D antenna array configurations have attracted much more research interest in array signal processing [22][23][24][25][26][27][28]. Most of those 3D arrays above are constructed from regular structure (i.e., cubic, cylinder), extending the planar array (i.e., URA, UCA) or configuring the virtual 3D array based on the planar array. Even though those special 3D arrays increase the elevation angle estimation accuracy, their array apertures are still relatively small since the physical size of static arrays are restricted. Moreover, conventional investigations in the DOA estimation require that the number of sensors is more than the number of receiving signals, which increase hardware cost and system complexity.
In contrast, in order to investigate the compromise between hardware cost and signal processing time, time-variant arrays whose element positions are changed over time are examined by time-divided sampling rather than simultaneously sampling as static array. Many researchers have reported using time-variant arrays to improve DOA estimation performance [29][30][31][32][33][34][35]. Instead of using a set of different elements to process the incident signals, the time-variant array can only use one or a small number of moving elements to construct virtual antenna arrays. Wan et al. proposed a method of combining the characteristics of arbitrary virtual baseline to construct virtual 3D array [29]. However, the number of sub-array elements is too small so that it cannot have high resolution. Liu examined a rotating long baseline interferometer whose length is much larger than one wavelength to estimate 2D DOAs by constructing virtual 2D circular arrays [30]. However, the 2D circular array has limited elevation aperture and still cannot well estimate the elevation angle.
The MUSB aperiodic array reconstructed from swarming UAVs proposed herein has aperture dimensions in both azimuth and elevation directions, which increase the accuracy of both azimuth and elevation angle estimation. Corner and Lamont proposed a parallel simulation of UAV swarm scenarios [36], and Saad et al. reported a testbed of vehicle swarm rapid prototyping [37]. Recently, many researchers investigated the methods and impact factors of designing the robust MUSB antenna arrays for signal collection platforms [38][39][40][41]. However, those works are limited in MUSB array constructing investigations including UAV positional precision, turbulence of the environment, micro-UAV swarm algorithm, and swarm-based real-time data collection. In this paper, 2D DOA estimation using the MUSB aperiodic array is provided, a mathematical model of the MUSB data collection system for signal processing is offered firstly and the impact of associated parameters on DOA estimation accuracy and convergence in this model are analyzed. The MUSB arrays have advantages of large aperture, large interelement spacing, no shadowing effect, low mutual coupling effect and large spatial sampling data from different locations in free space.
In practice, the MUSB system for DOA estimation requires low snapshots and might be applied in low SNR scenarios. However, the subspace-based techniques require adequate SNR and snapshots to guarantee good performance. We utilize the iterative method to lower the noise floor by multiplying the current MUSIC spectrum and previous spectrum for each iteration. The details of the iterative-MUSIC algorithm will be presented in Section 5.
This paper mainly contributes to the array signal processing area in three aspects. First, the mathematical model of the MUSB aperiodic array data collection system is introduced to demonstrate the fundamental DOA estimation impact factors and performance. Second, the CRB associated with DOAs in the presence of the gain and phase coefficient in the system is derived to reveal some direction-finding properties such as the global convergence, snapshots, SNR, and the number of arrays. Third, a successive DOA refinement procedure with iterative-MUSIC algorithm is provided based on the reconstructed arrays and spectrum to meet the requirement of high-precision DOA estimation. The rest of this paper mainly consists of six sections. Section 2 introduces the MUSB mathematical model; Section 3 derives the CRB associated with source DOAs in the presence of the sensor gain and phase coefficient; Section 4 proposes performance analysis of the MUSB system for the single-emitter case; Section 5 gives the algorithm applied in this paper; Section 6 provides simulation and measurement results in different scenarios and Section 7 concludes the paper.
Glossary of notations is listed below: C k×p = the space of k × p complex-valued matrices; E = expectation operator; A ij = the i, j element of a general matrix A ∈ C k×p ; A T = the transpose of A ∈ C k×p ; A H = the conjugate transpose of A ∈ C k×p ; Re (A) = the real part of A ∈ C k×p ; Im (A) = the image part of A ∈ C k×p ; tr (A) = the trace of A ∈ C k×k ; det (A) = the determinant of A ∈ C k×k ; A B = the Schur-Hadamard matrix product of A, B ∈ C k×p , defined by [A B] ij = A ij B ij ; A ⊗ B = the Kronecker matrix product of A, B ∈ C k×p , defined by z ∼ cN(µ(α), ζ(α)) = the complex Gaussian distribution of the complex random vector z with mean µ and variance ζ, and α is a real-valued parameter vector that completely and uniquely specifies the distribution of z (see [42]).

Swarming UAV Synthetic Aperture
A swarming UAV synthetic aperture was presented in our early published paper [43]. Figure 1 shows a graphical representation of a UAV swarm as it morphs in time (iteration I in this paper). Each of the M agents in the swarm has a location, orientation, and trajectory. Notionally, these have position P m,i (r, θ, φ), where m is the agent's index and i is the iteration. During swarming, the agents undergo rotations and translations, where a dual quaternion framework provides a convenient mechanism to handle this behavior. This motion rotates the agents' local (u, v, w) coordinate systems that describe the spatial orientation of their antenna radiation pattern with respect to (w.r.t.) the global coordinate system and the direction to the incoming signal of interest S(θ n , φ n ), which is the nth source. The collection of these measurements over iteration creates a synthetic aperture that can be used to calculate the parameters of interest (θ n , φ n ). Notionally, K is independent data sampled for each agent in each iteration and K is usually called a "snapshot".

Signal Model
Friedlander and Weiss presented a mutual coupling model in the presence of sensor mutual coupling, gain, and phase uncertainties [44]. We neglect the mutual coupling effect in the signal model since the spacing of aperiodic array reconstructed from swarming UAVs is expected to be much larger than a wavelength. Consider that an arbitrary array of M elements receive N uncorrelated incident signals in the far-field demonstrated in part 1. Thus, the received signal for the i-th iteration can be represented as where jω 0 τ i,mn and m = 1, 2, · · · , M; n = 1, 2, · · · , N; i = 1, 2, · · · , I. Then, A i . Note that the sensor gain g i,m and phase ψ i,m change w.r.t. element location based on orientation of UAV; τ i,mn changes w.r.t. location; S i is constant; W i may change w.r.t. location, velocity of UAV, and environment. Since the sources we consider here are in the far field from the observing array. It is easy to find that τ i,mn can be represented by τ i,mn = −d i,mn /c, and then where d i,mn is the distance from origin (reference sensor) of the coordinate to the m-th sensor in the direction of the nth source for the i-th iteration, c is the propagating velocity in free space, (x i,m , y i,m , z i,m ) are the coordinates of the m-th sensor for the i-th iteration, (θ n , φ n ) are the DOAs of the n-th source in the sphere coordinate. From Equations (5), the matrix A can be obtained by where λ is wavelength.

Swarming UAV Synthetic Aperture
In theory, for a static array, the steering vector A is considered invariant over different snapshots since the array geometry is invariant. Assume z ∼ cN(µ(α), ζ(α)), then, the m,n-th general formula of the Fisher information matrix (FIM) on the covariance matrix of any unbiased estimate of α is: where α m denotes the m-th component of α. The general formula has been presented in [45] and proved in [46]. Petre and Nehorai presented the deterministic and stochastic CRB in [46]. For deterministic CRB, the parameters, mean, and variance of the complex distribution are given by For stochastic CRB, the parameters, mean and variance of the complex distribution are given by Since the signals estimated cannot be known completely and the signals in the practical are stochastic, this paper considers the stochastic CRB.

CRB for the UAV Swarming System
Before deriving the CRB, we assume that both incident signals and noise are stationary and the ergodic complex Gaussian random process with zero mean and nonsingular covariance matrix is uncorrelated with each other. The columns of A = Γ A are linearly independent. An additional assumption is that the number of array elements reconstructed from swarming UAVs is greater than the number of sources. Therefore, the matrix of the steering vector has a full column rank.
The covariance matrices of signal, noise and observation vectors for i-th iteration are given by where A i is the steering vector for the i-th iteration and A i Γ i A i . It is useful to observe that if we let the sample data covariance matrixR i = 1 The log-likelihood function of K independent samples in i-th iteration of a zero-mean complex Gaussian random process X i (k) whose statistics depend on a parameter vector α is given by where Z denotes the constant term of the log-likelihood function, det(R) represents the determinant of the matrix R, and R i is the time-varying covariance matrix w.r.t. the iteration. Thus, the log-likelihood function of X(k) is: Therefore, the m,n-th elements of FIM are given by It follows that the FIM's submatrix F mn for the UAV swarming data collecting system can be obtained by summing the single-iteration F i,mn of FIM over the iterations. Furthermore, the FIM's submatrix F i,mn can be obtained by multiplying the single-snapshot F i,mn and the number of snapshots. Thus, we only need to know the single-snapshot F i,mn for the i-th iteration. The problem of the major interest is the estimation of the incident angles of the sources. Expression of CRB for 1D DOA of each iteration for the present problem is listed in [47]. The CRB of 2D DOA estimation with arbitrary array for the i-th iteration, which can be considered as an arbitrary static array, presented in this paper is given in the Appendix A.

Analysis of Single-Emitter Case
In this section, we investigate more details of the MUSB data collecting system with single-emitter case. The unknown parameters we consider here are the 2D DOAs (θ, φ). Assume the source variance is P, the noise variance is σ 2 , the snapshot for each iteration is K, and the iteration is I. From the Appendix A, we have the formula of FIM w.r.t. 2D DOAs in the i-th iteration of the MUSB system.
The mth element of the steering vector for the i-th iteration is given by where γ i,m is the gain and phase parameter and can be represented as γ i,m = g i,m e − jω 0 ψ i,m . Taking the derivative of a i (θ, φ) w.r.t. θ and φ for the i-th iteration, we obtain Therefore, it is straightforward to verify that Substituting Equations (14), (15), and (17), we can obtain Using the same derivative procedure, we can obtain Furthermore, Therefore, Assume where SNR = P/σ 2 . Thus, summing over iteration, we obtain the CRB where F θθ = 8Kπ 2 If we ignore the sensor orientation (i.e., Γ = 1) and let K = 1, then A H i A i = M and the FIM w.r.t. 1-D DOA is given by , which coincides with the results in [47].

UAV Parameters
Assume there are M swarming UAVs and each UAV swarms in a cylinder region (radius r, h = 2r).
Each UAV has an initial location (x, y, z) in the swarming region with a vector velocity → V i,m in each iteration and the initial locations of the M UAVs are considered as the first iteration. The swarming short distance in each iteration is represented as a vector → d i,m . The scalar quantity can be represented as d, which uses the following relation: where β t is a uniformly distributed random number matrix between zero and one, λ is wavelength and p is the coefficient determining the distribution mean of swarm distance. Figure 2 shows the distribution of short distance with mean µ = 0.25 wavelength. Here, the range of d also depends on the speed of UAV and data sampling interval and can be configured by the customer. We can obtain Figure 2. Distribution of short distance that one unmanned aerial vehicle (UAV) swarms.

Data Processing and Algorithm
When UAV swarms, a significant amount of data/information may be sampled due to the number of swarming UAVs and number of iterations. Herein, the location of each UAV is considered as one data point. Thus, when M UAVs morph I times, we have M * I data points which will be used to reconstruct virtual (e.g., synthetic) 3D aperiodic arrays to calculate the MUSIC spectrum and estimate DOAs. The number of data points (represented as N d ) is equivalent to the number of elements for a static array.
One challenge is how to choose N d for computing the MUSIC spectrum at each signal processing iteration (represented as p-iteration). We would like to use more N d in each p-iteration since more N d used for MUSIC spectrum calculation each time means that more array elements are used for data processing in a static antenna array, and are more accurate for DOA estimation. However, more data processing points cause higher calculation cost. Thus, it is necessary to compromise N d at each p-iteration and computational complexity.
Another problem is how to well choose data collecting and processing methods. Wide range of methods can be deployed to collect, concatenate, and process data collected from time-dependent measurements. MUSIC algorithm requires a minimum of three unique samples collected from these time-dependent measurements to provide DOA estimation. The required data samples can be obtained from agents in the swarm (e.g., one UAV locally collecting three samples at different times, or three spatially distributed UAVs collecting one sample simultaneously). We need to consider spatiotemporal distribution of samples (such as sampling rate w.r.t. wavelength, trajectory and velocity of UAVs, orientation of sensors, etc.) and iterative processing of measurements (use all data collected, truncate or applying a moving window, etc.). The iterative processing method with a moving window is called the iterative-MUSIC used in this paper. Figure 3 shows the data processing schematic for three data points in each iteration within the MUSB system. The left part in Figure 3 shows a UAV swarm and sampling method, and the right part shows the algorithm and data processing method. When UAV swarms to a certain location, we will sample K times and each snapshot takes time d t . After taking K snapshots, the program sets up a data point. When N d = 3, the program computes the stand MUSIC spectrum using three data points and stores the result. Then when UAV swarms again, we accumulate the current data point and two previous data points to calculate the MUSIC spectrum. Then we multiply the current and previous MUSIC spectrum at each p-iteration where we obtain the iterative-MUSIC spectrum to reduce noise level and improve DOA estimation performance. Note that if the p-iteration is too big, the value of spectral points will be very small and might be taken as zero. If so, DOA cannot be estimated and we may use dB instead of a number at that situation. Then, we refine and estimate DOAs and use predefined DOA estimation precision criteria to stop the process.

MUSIC Algorithm for MUSB Array
As stated in part 2 of this section, the number of reconstructed MUSB arrays is N d . Here, we rewrite the covariance of signals listed in Part 2 of Section 3 and give the first p-iteration (slightly different from UAV swarming iteration) of the iterative-MUSIC algorithm.
Rewriting the data model (1), we obtain where X{k} ∈ C N d ×1 are the vectors of sampled data, S{k} ∈ C N×1 are the source signals, Γ ∈ C N d ×N d are the gain and phase of sensors and A ∈ C N d ×N are the regular steering vectors. Thus, the covariance of X(k) is Define R s is N d × N d matrix with rank N. Therefore, it has N d − N repeated eigenvectors corresponding to the minimum eigenvalues σ 2 . Let e i be such an eigenvector so that R s e i = 0 or A H e i = 0, thus, N d − N eigenvectors e i corresponding to the minimal eigenvalues are orthogonal to each of N signal columns of A = Γ A, proved in [5]. N d − N dimensional subspace spanned by the noise eigenvectors is defined as noise subspace and N dimensional subspace spanned by incident signal mode vectors is defined as signal subspace. Let Q n be N d × (N d − N) matrix of noise eigenvectors, then the MUSIC spatial spectrum function is given by Then, the search spectrum peaks in the range of θ and φ, and the peak spectrum points we obtain are the estimation of arrival angles of incident waves.

Convergence Check
The algorithm performs the calculation until the system converges. The convergence can be guaranteed since the estimated DOA is a convergent series.
When the signal is covered by a high noise level, the estimated DOAs might be far from the ground truth and cannot be judged for convergence. However, as the iteration increases, the noise level is reduced and estimated DOAs are converged gradually. The Equation for judging convergence is given by where I m is the number of p-iteration of the iterative-MUSIC algorithm and ε is the preset threshold. The numerical simulation is given in part 1 of Section 6.

Computational Complexity Analysis
The orders of computational complexities of conventional spectral MUSIC [5], iterative-MUSIC, and eigenstructure-based algorithm with array interpolation (denoted by array interpolation) [48] are compared in Table 1. In this table, we assume the interpolated number of the array interpolation technique is equal to the number of actual sensors (M 1 = N d ) and the total number of angular sectors is denoted as I θ I φ . Those methods in Table 1 include the eigendecomposition step represented by the term O N 2 d N and the computation of J θ J φ samples of the MUSIC null-spectrum function represented by O J θ J φ (N d + 1) (N d − N) , where J θ and J φ stand for the search numbers along the directions of θ and φ. J θ J φ stands for the total search numbers for each sector of array interpolation algorithm. The iterative-MUSIC at each iteration may have very low snapshot (K = 1) different from the traditional spectral MUSIC algorithm which requires high snapshots. The computation complexity of the iterative-MUSIC algorithm is a little bit lower than the array interpolation at each iteration (J θ J φ ≈ I θ I φ J θ J φ >> N d ), but it is higher in the whole UAV swarming period. However, the array interpolation algorithm will be very complex when applied for 3D random arrays. Table 1. The orders of computational complexities of real-valued operations.

Simulation and Measurement Results
In this section, several groups of simulations will be carried out to demonstrate the performance of the presented distributed directional finding system in this paper. As the framework of the MUSB system established in this paper is mentioned for the first time, we focus mainly on analyzing the impact of various factors on the feasibility of the MUSB system. Practical measurements in the lab are also given in part 3 to show the DOA estimation performance in practice.
The wavelength of signal is fixed at 1 meter (m) and the simulation in each scenario is repeated 500 times. The elevation angle of the source emitter is 85 • and the azimuth angle is 270 • . As one UAV swarms till N d = 3, the iterative-MUSIC algorithm begins to search (0 • , 179 • ) space for elevation angle and (0 • , 359 • ) space for azimuth angle with 1 • interval to form the overcomplete MUSIC spectrum for each p-iteration. Then, UAV keeps swarming and the iterative-MUSIC algorithm keeps computing the MUSIC spectrum before the precision of DOA estimation is satisfied. When the preset threshold is satisfied, UAV stops swarming and the reconstructed process of phased arrays is terminated. The refined DOA estimations are obtained by scanning the reconstructed signal peaks from the iterative-MUSIC algorithm with 0.1 • step during the refinement procedure 10 times.
Moreover, the speed of UAV will influence the snapshot at each location where the system samples the source emitter. The snapshot at each location should be very low if the UAV swarms very fast and the snapshot can be high when the swarming rate of UAV is low. The speed of UAV in this paper will be represented by the distance between two iterations of swarming UAVs shown in Figure 2.
The joint root-mean-square error (RMSE) of incident signals is used for statistical DOA estimation precision evaluation, which is defined as where W is the number of Monte Carlo simulations, (θ w n , φ w n ) represent the actual DOAs of the nth signal, and θ w n ,φ w n represent the estimated DOAs of the nth signal in the w-th simulation. First, the system convergence will be studied in a typical scenario. Second, the DOA estimation performance using the iterative-MUSIC algorithm will be compared with CRB in various scenarios. Finally, DOA estimation performance in practice will be investigated.

DOA Estimation Performance
The performance of DOA estimation depends on multiple factors such as SNR, K, N d , velocity of UAV and number of iterations. Furthermore, the performance also depends on the distinctness of the array geometries due to the diversity of different observations at different time instants. In this section, we only consider the single-emitter case. When a single source is present, a typical scenario is set in which UAV swarming short distance mean µ = 2λ, the Monte Carlo simulation number and angles of the source are the same as before. Figure 5 depicts the results predicted by the stochastic CRB derived in Section 4. The joint RMSE of elevation and azimuth angles obtained from the iterative-MUSIC algorithm together with CRB from Section 4 are shown in Figures 6 and 7.
The extreme case is K = 1, N d = 3 (2D DOA estimation requires a minimum of three unique samples). The fixed settings and changed settings are listed as: Figure 5a Figure 5 give the lower bounds of the MUSB system. Scenarios (a)-(d) in Figure 6 show the performance of MUSB array without element rotation (sensor gain g = 1 and sensor phase ψ = 0 and scenarios (a)-(b) in Figure 7 show the impact factors with the element rotation of MUSB array (sensor gain and phase coefficients have certain deviations).
From Figures 6 and 7, we can find that Figure 6a shows the DOA estimation performance of the MUSB system will increase with increasing SNR; Figure 6b shows that the system can estimate DOAs even when snapshot K = 1; Figure 6c shows that when the UAV speed increases, the precision of DOA estimation increases; Figure 6d shows that the precision of DOA estimation increases with increasing elevation angles. Figure 7 shows the performance of the MUSB array with UAV rotating associated with sensor gain and phase varying. Figure 7a shows that the DOA estimation is not significantly impacted as the standard deviation of sensor gain increases; Figure 7b shows the precision of DOA estimation decreases significantly when the standard deviation of the sensor phase increases to a certain value.
In Figure 7b, we find that the stochastic CRB is flat for the single-emitter case as the phase standard phase varies because Equation 16 shown in Section 4 cancels the sensor phase errors by multiplying the steering vector A and the conjugate transpose of A, while the iterative-MUSIC algorithm does not have the advantages CRB has utilized. Equation 39 shows the advantages of CRB for the single-emitter case and Equations 16-21 give the derivation process with advantages of sensor phase error cancel.
where ϕ i,m includes the sensor phase deviation and phase difference from sensor m to relative phase center in the i-th iteration.

Measurement
In the experiment, the test fixture provides a convenient platform to study this morphing in time (using sixteen elements). Randomly positioned rectangular patch antennas designed for 2.45 GHz are used with a randomly morphing volume provided by a moving platform named as "Medusa" to estimate DOAs. Figure 8 shows that one source in the far-field transmits and receiving sensors receive the signal with the vector network analyzer (VNA) measuring the S 21 of transmitting and receiving antennas. Figure 9 shows the measured MUSIC spectrum with 16 spatially distributed elements as a static volumetric random array in Medusa with element rotation (Rotate 0-45 degrees around x, y, z axis randomly, no UAV swarming). Figure 10 shows the measured azimuth and elevation errors gradually converge from around 10-20 degrees to around 1-2 degrees as the iteration increases with UAV swarming. In Figure 10, we take the extreme case (N d = 3, K = 1).

Conclusions
This paper establishes a MUSB data collection framework for the first time, which makes it possible to realize source estimation with low snapshot under low SNR environment in a UAV swarming period. Theoretical and experiment results are given to reveal the performance of the MUSB phased array system used for 2D DOA estimation, which supports the feasibility of the system. The iterative-MUSIC algorithm is applied for the framework and it can estimate the DOAs efficiently only with one snapshot in each iteration when UAV swarms very fast. The UAV speed controls the structure of the reconstructed aperiodic phased arrays from the MUSB system and the DOA estimation precision is increased when the distance between the two iterations of swarming UAV is increased. The impact of known sensor gain errors and phase errors from the UAV rotation for the DOA estimation performance are also investigated. Practical experiment results match the theoretical expectation of the MUSB system using the iterative-MUSIC algorithm. Our results will benefit future research on performance analysis and optimal design of time-varying antenna arrays based on the UAV swarm. It is also interesting to extend the results when position errors are present in the future.

Conflicts of Interest:
The authors declare no conflict of interest.