A new model for DOA estimation and its solution by multi-target intermittent particle swarm optimization

Currently, the widely used methods for direction of arrival (DOA) estimation were constructed based on the subspace, such as Multiple Signal Classification (MUSIC) and Estimating Signal Parameter via Rotational Invariance Techniques (ESPRIT), which required that the number of sources is known beforehand. In this paper, a new method based on the Vector Error Model (VEM) for estimating the DOAs was proposed, which do not need the sources number in advance. The comparison of the performance between the VEM and the MUSIC model for DOA problem was given to demonstrate the effectiveness of our method. The algorithm of multi-target intermittent particle swarm optimization (MIPSO) was adopted to solve the VEM, and the performance of the VEM-MIPSO method was analysed through simulations for a uniform linear array and an L-shaped array respectively. The results showed that: (1) the VEM was an effective model to solve the DOA estimation without prior knowledge of the sources number; (2) the MIPSO was an efficient algorithm to solve the DOA estimation with high precision.


Introduction
In the past couple of decades, the problem of direction of arrival (DOA) estimation for narrowband sources has been studied extensively because of its wide application in many fields such as radar, navigation and wireless communication (Lee et al., 2018). The subspace methods of multiple signals classification (MUSIC, proposed in Schmidt and Schmidt (1986)) and Estimating signal Parameter via Rotational Invariance Techniques (ESPRIT, proposed in Yuen and Friedlander (1996)), are the most prominent approaches for DOA estimation with high resolution and simple structure. The theory of the MUSIC methods is based on the assumption of partitioning the observation space into two orthogonal subspaces of signal subspace and noise subspace. Actually, these methods cast the DOA estimation as an optimization problem objective function which involves the empirical estimate of the signal subspace projection matrix (Zhou, Huang, et al., 2018). However, the number of sources was assumed to be known or estimated by an information theoretic criterion, which could not be satisfied usually. Some source enumeration techniques, such as the Akaike information criterion (AIC), the minimum description length (MDL), have been proposed to solve the problem faced by the approaches based on MUSIC, but the performance CONTACT Peichao Zhao 781961774@qq.com degrades heavily if the signal to noise ratios (SNR) is low or the number of samples (snapshots) is small (Djuric, 1996). In another aspect, researches for quantifying performance loss of MUSIC-based DOA estimation has been addressed, which turn out to be a function of the number (M) of sensors, number of snapshots (N) and signalto-noise ratio (SNR), respectively. However, the statistical performance is only characterized in the situation where N tends to infinity while M remains finite. When M is probably close to or even larger than N, strong discrepancies between the sample eigenvectors and their population counterparts will occur (Thiergart & Habets, 2013). Therefore, a method which could finish DOA estimation without knowing the number of sources is needed.
Some DOA estimation algorithms without source number have been proposed in Qi et al. (2004), Jin et al. (2008, and Zhang and Ng (2010). The method proposed in Qi et al. (2004) works with a threshold for pre-estimation of the number of sources, which was not addressed. An extension of the reversible jump Markov chain Monte Carlo method is mentioned in Jin et al. (2008), which is very inefficient and computationally intensive. A MUSIC-like method was introduced in Zhang and Ng (2010), which was constructed based on the framework of beamforming.
In this paper, a method based on the Vector Error Model (VEM) was adopted to solve the problem of DOA estimation, and an algorithm named multi-target intermittent particle swarm optimization (MIPSO) was used to solve the model based on the VEM. The rest of this paper is organized as follows: Section 2 introduces the VEM for the DOA estimation problem and the solving algorithm of MIPSO. Several simulations were carried out in Section 3 to demonstrate the performance of the VEM-MIPSO method for DOA estimation. Finally, conclusions were drawn in Section 4.

Solving method
In this section, the MUSIC model for the DOA estimation was mentioned firstly. Then, the new algorithm based on the VEM was proposed. (Zhou, Hu, et al., 2018) Consider the situation as shown in Figure 1 (a), where the ith source impinging from the far-field at an angle of θ i onto a Uniform Linear Array (ULA) of M identical omnidirectional equidistant sensors with an inter-sensor space of d. The received signal at the mth sensor is a superposition of time-shifted versions of all the sources corrupted with noise as following

Music model for DOA estimation
where c 0 is the speed of propagation, K is the number of the sources, t is the sampling point for the time. Rewriting (1) in vector format we get where T with the superscript (·) T denoting the transpose is a complex matrix required by the sensors; S(t) = [s 1 (t), s 2 (t), . . . s K (t)] T with the sources number of K contains all the narrow-band far-field source signals; ; λ is the wavelength of the source; and N(t) = [n 1 (t), n 2 (t), . . . , n M (t)] T is an unknown noise matrix.
The aim of DOA estimation is to calculate all theθ i , i = 1, 2, . . . ,K based on (2), where the top mark of ∧ represents the estimated value. Generally, as a fixed sampling grid that uniformly covers the DOA range [−π/2, π/2]. An L-shaped array, as shown in Figure 1 (b), consists of two orthogonal (M + N)-element ULAs (denoted by x-ULA and y-ULA, respectively) in the xy plane, the narrow-band far-filed uncorrelated signals impinge on the array from the distinct directions with elevation and azimuth angles {(θ i , φ i )}. Different from (2), the steering matrix for the Lshaped array has the form shown in (3).

DOA estimation model based on VEM
In order to avoid the calculation for the complex, we redefine (2) as follow is the real part of the complex signal matrix from M sensors' array without the superscript of (·) T ; real(·) represents the real part of a complex value; T is the transpose of the real part of the steering matrix A(θ ); and V(t) is the transpose of the real part of the noise matrix N(t).
According to (4), the task of DOA estimation become a problem to judge whether a curve c(θ) with the parameter ofθ exists in the matrix C(θ ) as a row vector only based on the measured data D(t). In order to obtain this c(θ), a matrix W is introduced so that . . .
(5) where the note of → represents that the Y approximates the C(θ ). Therefore the objective function is given as where || · || 2 2 is the 2-norm of a vector; ∀i means for all parameters. Solving (6), a model named Vector Generated Model (VGM) described in (7) is obtained, which could generate a vector y T i according to the known vector c T (θ i ). Equation (6) is called VEM, which calculated out the error between the known vector c T (θ i ) and the generated vector y T i .
whereD * is a matrix generated from the matrix D. Please refer to references , , Cui et al. (2015) or Appendix 1 to find the calculation process from (6) to (7). According to the definition of VEM in (6), if the known vector c T (θ i ) exists in the matrix C(θ ), the generated vector y T i will equal to the known vec- where real(·) represents the real part of a complex value. Similarly, we define the VEM shown in (9) for the L-shaped array.

Solving algorithm
First, let us consider the problem described in (6), then (10) will be solved similarly. There are many algorithms to find the optimal solution for (6). However, there should be more than one solution for (6) needing to be obtained, and the number of the solution is unknown which corresponds to the number of the signal Here the multi-target intermittent particle swarm optimization (MIPSO) algorithm is adopted to calculate the optimal parameters where θ i is the current position (value) of the ith particle (parameter), θ + i is the next position of the ith particle, θ i is the current velocity (increment) of the ith particle, θ + i is the next velocity of the ith particle, b i is the best position of the ith particle, t i is a particle around the ith particle within a small scope δ, and ε(t i ) is the smallest value in this scope δ. ω is the inertia coefficient, c 1 and c 2 are two accelerated factors which are always set as 2, r 1 and r 2 are two random constants whose value are between 0 and 1. Here, t i replaces the global best particle in the traditional PSO algorithm. During the calculation, every particle will find it's t i and will approach this t i dynamically. If the number of these {t i } L i=1 and their fitness ε(t i ) do not change during a certain step, or if the maximum step is reached, the calculation iteration will stop. In order to keep the ability of exploitation for all the particles to find their t i , these t i are searched intermittently (every 20 steps). After the iteration stops, all the particles will converge around several t i . The final solution where N ≤ L, according to two criteria: (1) whether the ε(t i ) is smaller than a threshold and (2) whether the values of the t i locate within the definition domain. Please refer to reference (Cui et al., 2015) for detailed information.
For objective function for the L-shaped array shown in (10), the MIPSO algorithm is designed as where the parameter of θ i is replaced by (θ i , φ i ). Figure 2 gives the process of the calculation based on VEM and MIPSO. In order to explain it simply, we use the parameter θ i for the ULA in Figure 1(a). Please replace the parameter of θ i by the parameter of (θ i , φ i ) for the L-shaped array.

The VEM-MIPSO method
Firstly, the angle parameters θ i were initialized randomly with the number of P. Then, reference curves c T (θ i ) were generated according to the θ i . Input the c T (θ i ) to VGM, y T i were calculated. The errors ε(θ i ) were generated according to c T (θ i ) and y T i . Finally, the parameters θ i were adjusted based on the errors ε(θ i ) until all the parameters θ i converged around several ones {t i } L i=1 . The criteria mentioned above is used here to pick the final solutions where K is the number of the sources, and much smaller than P.
For calculation convenience, the errors between y T i and c T (θ i ) are defined as where the subscript of P is the number of the parameters calculated by the MIPSO algorithm simultaneously. According to the definition in formula (14), the ε(θ i ) will decrease when the c T (θ i ) and y T i approach the same shape.

Simulation and analysis
In this section, simulations to investigate the performance of the VEM-MIPSO method proposed in this paper are conducted. The program for the DOA estimation was written by Matlab R2014b, which could be obtained through the email of clzh0308@126.com. We consider ULAs with 50 sensors, and inter-sensor spacing of d = λ/2. Three echo signals with DOA of θ 1 = 10 • , θ 2 = 30 • and θ 3 = 60 • , respectively, were generated by program. For the Lshaped array, 50 sensors are used for both x-and y-axis respectively. The pairs of elevation angle and azimuth angle are set as (θ 1 , φ 1 ) = (10 • , 15 • ), (θ 2 , φ 2 ) = (30 • , 25 • ) and (θ 3 , φ 3 ) = (60 • , 35 • ). The snapshot for both arrays was set as 10. All the experiments are carried out in Matlab R2014b on a PC with Intel E33-1231 v3 CPU and 16GB of RAM.

Analysis of the model
The performance of the VEM model for DOA estimation was firstly analysed in this section. Figure 3 gives the comparison of the performance between the VEM model and MUSIC method under different level of noise jamming for ULA. Table 1 lists the errors calculated by VEM model under different level of SNR. From Figure 3 and Table 1, the following comments could be given.
(1) Similar to MUSIC method, the VEM could give featured values at the direction of arrival. The difference is that the VEM gives local minimum errors at the target angles, while the MUSIC method gives local maximum values in the spatial spectra at these angles.
(2) The accuracy of the VEM is higher when the noise level becomes severe. In Figure 3(b), there is a deviation at the angle of 60 o when the SNR equals to 1 dB, while the VEM gives all the local minimum value at correct angles. Particularly, calculation error occurs   Figure 3 gives the comparison of the performance between the VEM and MUSIC method under different level of noise jamming for L-shaped array. Table 2 lists the errors calculated by VEM model under different level of SNR. From Figure 2 and Table 1, the following comments could be given.
(1) Both of the two methods could give featured values at the direction of arrival. The difference is that the VEM gives local minimum errors at the solutions, while the MUSIC method gives local maximum values in the spatial spectra at these angles.
(2) The resolution of the MUSIC method for L-shaped array is higher than that of the VEM model. From the Figure 4 (b), we can see that only the target coordinates have a local maximum value, and the others points own little values without severe fluctuation. While the errors calculated by VEM model, as shown in Figure 4 (a), have a bigger fluctuation among the points.
(3) The model based on the VEM has quicker speed (0.06 s), to draw the whole error map than the MUSIC method (0.1 s). (4) From Table 2, we can see the minimum values at certain angles remain small level under 0.3. Therefore, the value of 0.3 was selected as the threshold during the calculation of MIPSO algorithm.

Analysis of calculated results
The performance was evaluated by the root mean square error (RMSE), which is calculated by formula (16). For the ULA, the comparison is carried out between VEM-MIPSO method and the root-MUSIC method. The RMSE comparisons for ULA under different SNR, sensors' number and snapshots are list in Figure 5 (a-c), respectively.
(1) As shown in Figure 5(a), the experiment for both the two methods under different SNR situation is executed while setting the sensors' number as 50 and the snapshot as 10. Both of the errors generated by the VEM-MIPSO method and the root-MUSIC method could remain at a low level, even when the noise existing in the data is severe. But the time used by the VEM-MIPSO (with average of 0.98 s) is much shorter than that used by the root-MUSIC method (with average of 2.2 s). (2) Figure 5(b) gives the experimental results when increasing the sensors' number from 20 to 60 by the step of 2 under the fixing snapshots of 10 and SNR as 20. From the trend of the RMSE, the accuracy from (1) The RMSE comparison for L-shaped array with different noise levels when setting the sensors' number as 50 and the snapshot as 10 is given in Figure 6, from which we can see that although the time used for the VEM-MIPSO method is longer (with average of 1.2 s) than that for the MEMP method (with average of 0.45 s), the RMSE of VEM-MIPSO model is much lesser.
(2) In Figure 7, we changed the sensors' number of the array from 40 to 60 by step of 1 while fixing the snapshot as 10 and the SNR as 20. Different from the situation for ULA in Figure 5(b), the sensors' number adopted in the array has little effect for the MEMP method, while the RMSE for the VEM-MIPSO method become unstable when the sensors' number decreases down to a certain value. (3) Finally, the snapshots of the data were adjusted while fixing the sensors' number as 50 and SNR as 20. The RMSE comparison is described in Figure 8. Similar to the situation for ULA in Figure 5(c), the snapshots  involved in the calculation has little effect for the MEMP method, while the RMSE for the method based on VEM increases significantly when the snapshot goes up to a certain value under a fixed sensors' number. So it is better to set a suitable ratio between the sensors' number and the snapshot before using the VEM-MIPSO method.

Conclusion
This paper proposed a new method to estimate the DOAs based on the VEM, which had a totally different view for the previous ones. An algorithm named MIPSO was adopted to solve the model without knowing the number of the signal sources. And several simulations with deep analysis were given to illustrate the efficiency and advantages of the proposed method. Following conclusions were drawn:

Efficiencies
The method based on the VEM maps the solution domain of the DOA to an error space, which is totally different from the current methods mapping the solutions to a spatial spectra space. In this error space, only several solution parameters have the local minimum values which are much lesser than the values around them. This feature of the error space provided the possibility that the minimum value could be calculated by warm intelligent optimization algorithm. Considering the fact that multitarget (multi-angle) need to be found, an algorithm of MIPSO was adopted to find all the targets simultaneously.

Advantages
For the ULA, the method based on VEM and MIPSO algorithm could find the solutions with similar accuracy as the root-MUSIC method but much lesser time. For the L-shaped array, the proposed method could calculate the solutions with much smaller RMSE than the MEMP method.

Limitation and further work
The ratio between the sensors' number and the snapshot should be set at a suitable value to guarantee the efficiency of the VEM-MIPSO method. According to the simulations in this paper, the good accuracy could be obtained with sensors' number of 47 and the snapshot of 10.
For the L-shaped array, the VEM-MIPSO method need more time to finish the calculation than the MEMP method, which was caused by the process of the MIPSO. This limitation could be solved by improving the MIPSO algorithm to accelerate the speed of the convergence of the particles.
In this paper, the situations of the ULA and the Lshaped array are analysed. And only the far-field narrowband signals are concerned. Further experiments could be tried for different types of sensors' array and signals.
The resolution of the VEM-MIPSO method among different DOAs could be refined by adjusting the parameter δ in the future work.

Disclosure statement
No potential conflict of interest was reported by the author(s).