Signal Subspace Reconstruction for DOA Detection Using Quantum-Behaved Particle Swarm Optimization

: Spatial spectrum estimation, also known as direction of arrival (DOA) detection, is a popular issue in many ﬁelds, including remote sensing, radar, communication, sonar, seismic explo-ration, radio astronomy, and biomedical engineering. MUltiple SIgnal Classiﬁcation (MUSIC) and Estimation Signal Parameter via Rotational Invariance Technique (ESPRIT), which are well-known for their high-resolution capability for detecting DOA, are two examples of an eigen-subspace algorithm. However, missed detection and estimation accuracy reduction often occur due to the low signal-to-noise ratio (SNR) and snapshot deﬁciency (small time-domain samples of the observed signal), especially for sources with different SNRs. To avoid the above problems, in this study, we develop a DOA detection approach through signal subspace reconstruction using Quantum-Behaved Particle Swarm Optimization (QPSO). In the developed scheme, according to received data, a noise subspace is established through performing an eigen-decomposition operation on a sampling covariance matrix. Then, a collection of angles randomly selected from the observation space are used to build a potential signal subspace on the basis of the steering matrix of the array. Afterwards, making use of the fact that the signal space is orthogonal to the noise subspace, a cost function, which contains the desired DOA information, is designed. Thus, the problem of capturing the DOA information can be transformed into the optimization of the already constructed cost function. In this respect, the DOA ﬁnding of multiple signal sources—that is, the multi-objective optimization problem—can be regarded as a single objective optimization problem, which can effectively reduce the probability of missed detection of the signals. Subsequently, the QPSO is employed to determine an optimal signal subspace by minimizing the orthogonality error so as to obtain the DOA. Ultimately, the performance of DOA detection is improved. An explicit analysis and derivation of the developed scheme are provided. The results of computer simulation show that the proposed scheme has superior estimation performance when detecting signals with very different SNR levels and small snapshots.

The problem of DOA detection has been addressed by a large number of methods, as seen in the literature. Among them, high-resolution eigen-subspace algorithms are the most studied, and emerging algorithms are mostly based on the concept of the eigensubspace algorithm. Specifically, MUltiple SIgnal Classification (MUSIC) [12][13][14][15][16] and the Estimation Signal Parameter via Rotational Invariance Technique (ESPRIT) [17][18][19][20][21] are the most representative methods encountered in DOA detection.
The MUSIC method, which has a high resolution for signal source detection in ideal conditions, is applicable to any array configuration. The MUSIC algorithm needs enough snapshots to construct a sampling covariance matrix, which is used to replace the data covariance matrix, and DOA information can be obtained through a spectral peak search. However, due to small signal-to-noise ratios (SNR), lack of snapshots, and spectral peak diffusion [22], the performance of MUSIC deteriorates or becomes invalid, leading to missed detection and low detection accuracy. In addition, the spectral search process of MUSIC is essentially a multi-objective optimization problem. Thus, its performance may significantly decrease in a scenario involving signals with very different SNR levels, especially for low SNR thresholds and small snapshots.
The ESPRIT method has been considered a visible technique since its introduction due to its search-free estimation. ESPRIT exploits the subarray partition and rotational invariant of different subarrays to obtain high-resolution DOA information. Nevertheless, it suffers from a loss of array aperture to a certain extent, which results in performance degradation for DOA detection, especially for low SNR thresholds. In addition, the ESPRIT approach is based on a uniform linear array (ULA). Some additional preprocessing operation is indispensable when applying it to arbitrary antenna configurations, such as array transformation [23] or special array geometry's high-order cumulant calculation [24].
Another kind of DOA detection algorithm is the subspace fitting algorithm [25][26][27][28][29]; a typical one is the weighted subspace fitting (WSF) algorithm [30,31]. Its fundamental principle is that there is a fitting relationship between the subspace of the received data and the array manifold matrix, based on which an estimation of DOA information is obtained through solving a maximum or minimum likelihood function. When some prior information is known, such as the number of signal sources, the WSF algorithm is more adaptable to a low SNR and fewer snapshots compared with the subspace-class algorithm [32,33]. In addition, the WSF algorithm can effectively estimate DOA for both coherent and incoherent signal sources, while the traditional subspace-class algorithm is only valid for incoherent sources. Nonetheless, due to the nonlinear property of the likelihood function, the WSF algorithm needs a multi-dimensional joint search for multiple angles, which brings a high and unacceptable level of complexity to the computation, rendering it unsuitable for a practical engineering application.
From the above analysis, the subspace-class algorithm based on covariance matrix decomposition has poor adaptability to a low SNR and fewer snapshots, and the high computational complexity of a subspace fitting algorithm is often unacceptable in practical application. Therefore, it is necessary to conduct further research in DOA detection when the SNR threshold is low and the number of snapshots is insufficient.
It is known that in the case of a low SNR threshold, improving the performance of DOA detection is very difficult. The reason is that in real-world situations, the environment can be very noisy (e.g., as encountered in a battlefield remote sensing scout system [34]). To further improve the DOA detection performance when the value of the SNR and the number of spatial samples are small, in this paper, an augmented DOA detection method with signal subspace reconstruction using Quantum-Behaved Particle Swarm Optimization (QPSO) [35][36][37] is proposed, which has strong adaptability to a low SNR and a small sampling number. Firstly, we use less snapshot data to obtain a sampling covariance matrix under a low SNR, which is exploited to build a noise subspace through eigendecomposition. Then, we randomly select several directions from the target space to build a potential signal subspace with a signal steering matrix. A multi-dimensional cost function including the DOA information of interest is constructed with the aid of the orthogonality between the signal subspace and the noise subspace, through which a multi-objective optimization problem about the DOA finding of multiple signal sources is regarded as a single-objective optimization, such that DOAs of multiple signal sources can be included in a single spectral peak. Subsequently, the QPSO is considered in order to optimize the directions (cost function) to the modified steering matrix and signal subspace. After carrying out the above process, the performance of DOA detection becomes more robust. The detailed formula derivation and simulation experiments of the proposed scheme are provided. To the best of our knowledge, the method proposed in this paper has not been studied elsewhere in the literature.
The remainder of this study is arranged as follows. An array signal model is formulated and a subspace-based algorithm is introduced briefly in Section 2. In Section 3, a novel DOA detection approach, using the QPSO algorithm with signal subspace reconstruction, is proposed and related formulas are derived. In Section 4, simulation results are presented. Section 5 concludes the study.

Array Signal Model
The transmission process of a signal in a wireless channel is very complex. The array antenna is not always omnidirectional, such as in a cylindrical conformal array with directional antenna elements [38]. Sensor array error is also a non-negligible factor and includes antenna element pattern error, channel uncertainty, array element mutual coupling and sensor position error [39][40][41]. If these errors are not compensated for effectively, the performance of a high-resolution spatial spectrum estimation algorithm will deteriorate significantly or even fail. In addition, for different signal types (e.g., far-field signal and near-field signal, broadband signal and narrowband signal, colored noise and Gaussian white noise), the mathematical models of the array and signal processing steps will also display differences [42][43][44]. In this case, a practical parameter model is adopted after some consideration and for the purpose of simplicity. The following assumptions are about antenna arrays and received signals, which are applied to the following sections: (1) An array antenna is omnidirectional-that is, the elements in the antenna radiate uniformly in all directions of space. The complex gain of each element in the antenna array is the same without considering mutual coupling.
(2) The noise contained in the received signal conforms to Gaussian distribution with zero mean, and it is independent of the observed signal.
(3) The array receives plane waves from the far field of the signal source.
(4) All signals are narrow-band. Consider a ULA consisting of M (m = 1, 2, · · · , M) sensor elements on which P (P < M) incoherent signals s p (t) (p = 1, 2, · · · , P) are impinging with different angles θ 1 , θ 2 , · · · , θ p , · · · , θ P , where θ P is measured with respect to the y-axis. The model of the array signal is illustrated in Figure 1. Let the sensor at the coordinate origin be the reference and d be the spacing of adjacent sensor elements, and the phase difference between two adjacent sensors is defined as where λ (λ ≥ 2d) denotes the wavelength. The direction vector of such ULA at angle θ p can be expressed in the form where the superscript [•] T represents the transpose operation. Let A be an M × P array manifold matrix obtained by stacking the direction vectors a θ p of the ULA, i.e., The array output at the tth snapshot is computed as follows: where y(t) is the M × 1 measured data, s(t) = [s 1 (t), s 2 (t), · · · , s P (t)] T is a P × 1 signal vector, and the M × 1 vector n(t) = [n 1 (t), n 2 (t), · · · , n M (t)] T represents a Gaussian- distributed noise vector that is assumed to be spatially white. The covariance matrix of the y(t) is given in the following form: where H, σ 2 n , and I represent the complex conjugate transpose, the variance of the additive white Gaussian noise, and an identity matrix, respectively.
t n t n t n t n represents a Gaussian-distributed noise vector that is assumed to be spatially white. The covariance matrix of the   t y is given in the following form: where H ,  2 n , and Ι represent the complex conjugate transpose, the variance of the additive white Gaussian noise, and an identity matrix, respectively.
The modeling scheme of measured data.

Introduction of MUSIC Algorithm
The basic principle of MUSIC is to extract the signal subspace and noise subspace from the sampling covariance matrix of array output data, and then exploit the orthogonality of these two subspaces to estimate the signal parameters. The eigenvalue decomposition of data covariance matrix R is given by where S U is the signal subspace composed of eigenvectors corresponding to P large eigenvalues. N U is the noise subspace, which is expanded by eigenvectors corresponding to the remaining  M P small eigenvalues. S Σ and N Σ are diagonal matrices with signal and noise eigenvalues. Considering that the length of the receiving data vector is finite, the maximum likelihood estimation (MLE) of R is computed as where L is the number of snapshots. Thus, the spatial spectrum of MUSIC is derived in the following form by or Figure 1. The modeling scheme of measured data.

Introduction of MUSIC Algorithm
The basic principle of MUSIC is to extract the signal subspace and noise subspace from the sampling covariance matrix of array output data, and then exploit the orthogonality of these two subspaces to estimate the signal parameters. The eigenvalue decomposition of data covariance matrix R is given by where U S is the signal subspace composed of eigenvectors corresponding to P large eigenvalues. U N is the noise subspace, which is expanded by eigenvectors corresponding to the remaining M − P small eigenvalues. Σ S and Σ N are diagonal matrices with signal and noise eigenvalues. Considering that the length of the receiving data vector is finite, the maximum likelihood estimation (MLE) of R is computed aŝ where L is the number of snapshots. Thus, the spatial spectrum of MUSIC is derived in the following form by It is well known that the spatial spectrum will present a peak if θ equals the actual DOA [12]. Typically, the spectral peaks are searched within a certain angle range (e.g., θ ∈ [θ l , θ r ]), and the DOAs to be estimated are determined by the P largest spectral peaks. In other words, DOA information of different signal sources is included in different peaks, and the estimation process is performed in a serial manner. This problem is essentially a multi-objective optimization problem. Although the estimation accuracy of the MUSIC ap-Remote Sens. 2021, 13, 2560 5 of 12 proach is close to the Cramer-Rao Bound (CRB) [45][46][47] in the ideal case, the performance deteriorates rapidly in the case of a low SNR and few samples, especially for signals with very different SNR levels. Since the angles corresponding to the top P spectrum peaks are regarded as the true DOAs in the MUSIC algorithm, false alarm, missed detection, and precision degradation occur easily in the event of a low SNR and snapshot deficiency. Moreover, the MUSIC algorithm suffers from an excessive computational burden in multidimensional parameter detection due to its search-based strategy.

DOA Detection Approach Using QPSO through Signal Subspace Reconstruction
In this section, a direction-finding method through the reconstruction of the signal subspace with the QPSO algorithm will be described in detail. This method, as we will demonstrate, has excellent estimation performance under scenarios with a low SNR and snapshot deficiency.
In view of the above analysis, the constructed cost function leads to a possible multimodality problem, which cannot be solved efficiently through traditional mathematical methods. However, QPSO is a promising alternative. In a classical Particle Swarm Optimization (PSO) system [48][49][50], the actual search space is unable to cover the whole solution space because the speed of the particles is always limited, and the convergence of the particle motion relies on its trajectory. In addition, the PSO algorithm cannot guarantee global coverage with a probability of one, which is a drawback of this algorithm. In 2004, in line with the philosophy of PSO, Sun proposed the QPSO algorithm [35][36][37]. In Newtonian mechanics, the motion of each particle follows its own trajectory. On the contrary, according to the interpretation of quantum mechanics, the trajectory is meaningless in the consideration of the uncertainty principle [51]. Therefore, with a certain probability, particles can appear anywhere in the solution space, which indicates that the QPSO algorithm is able to converge to a global optimum value. Obviously, the QPSO algorithm is a probability-based PSO algorithm, and its key idea is to establish a potential field to constrain particles.
In an observation space, P different angles are randomly selected to form an initial vector Θ = θ 1 , θ 2 , · · · , θ p , · · · , θ P . Based on the signal subspace principle, we reconstruct a signal subspace, which is shown as A Θ = a θ 1 , a θ 2 , · · · , a θ p , · · · , a θ P (11) where A Θ is the array manifold matrix, and a θ p is the steering vector of θ p . To obtain a sound signal subspace, we build and optimize such a cost function due to the orthogonality between the two subspaces in the following: where · 2 denotes l 2 norm. As mentioned above, the cost function is a possible multimodality problem, and the QPSO will be employed for its good astringency and fast convergence speed to obtain the true DOA information of signal sources. Thus, Equation (13) is modified as In P dimensional space with a population size of N, the position vector of the nth (n = 1, 2, · · · , N) particle is Θ n = (θ n1 , θ n2 , · · · , θ nP ). The fitness value is calculated by substituting Θ n into the cost function min f Θ . The personal best position of the nth particle is , g ∈ {1, 2, · · · , N}, where g denotes the subscript of the particle in the global best angle. The mean of the personal best positions of all particles is where t denotes the current iterative times. In terms of the convergence behavior of particles in [52], every particle should converge to its local attractor Θ la = θ la 1 , θ la 2 , · · · , θ la N , which is Θ la Applying the Monte Carlo method, the particle position is determined by where α is called the contraction-expansion coefficient. Normally, the value of α is calculated as [37] where t max is the maximum iteration time. In general, w 1 = 0.5, w 2 = 1.
In summary, the basic steps are: (1) Initialize particle position vector Θ n and the best previous position Θ pbest n of each particle. (2) Calculate the mean best position mbest(t) using Equation (15).
(3) Determine the current fitness value by substituting Θ n into the cost function min f Θ and compare it with the particle's previous best value. If the current fitness value is smaller than the previous best value, replace the previous best value with the current value.   [53], in which only one parameter, namely the position vector, is used to describe the particle state. Therefore, the algorithm has a fast convergence rate, has few parameters to set, and is easy to implement.
It should be noted that although Θ opt n = θ opt n1 , θ opt n2 , · · · , θ opt nP is the so-called optimal solution of the particle "position" vector, θ opt n1 , θ opt n2 , · · · , θ opt nP is actually the desired P DOAs. Therefore, the proposed method transforms the multi-objective optimization problem in the traditional search-based DOA detection method, into a single-objective optimization problem, which can obtain the DOAs of multiple signal sources in a single spectral peak, avoiding false alarm and missed detection, and it performs well in the event of a few snapshots and a low SNR.

Experimental Studies
To compare the root mean square errors (RMSE) of the proposed approach with MUSIC and ESPRIT, simulation experiments were conducted in MATLAB (MathWorks.Inc. Natick, Massachusetts, USA). In this paper, the RMSE [54] criterion was defined as where P denotes the number of signal sources, N e denotes the number of independent experiments, θ p is the pth DOA of the incident signal source, andθ p (n) is an estimate of θ p in the nth Monte Carlo trial.
In the first experiment, the algorithm parameters were set as follows: there were 10 elements in the ULA, and the spacing between adjacent elements was d = λ/2. We assumed that three incoherent signals coming from 5 • , 10 • , and 5 • impinged on the ULA. The SNRs of the three signals were set to be −15 dB, and the snapshot varied from 10 to 50 with 10 intervals; for the MUSIC algorithm, we set the search step as 0.1 • . A total of 100 experiments were completed. Figure 3a shows the performance of the RMSE versus snapshots (signals with the same SNR). It was apparent that the proposed method outperformed the other methods. The RMSE of the proposed method became smaller as the number of snapshots increased. The principles of MUSIC and ESPRIT determine that their performance is highly dependent on the covariance matrix and the signal/noise subspace. Under the condition of a low SNR and small number of snapshots, their performance deteriorates significantly or even fails due to the poor accuracy of the covariance matrix of the array output and signal/noise subspaces. In addition, in Figure 3b, we also plot the cost function and the RMSE versus the number of iterations. The number of snapshots was fixed at 30. It was noticeable that the cost function and the RMSE decreased as the iterative process proceeded, and the method could converge within 45 iterations.

Experimental Studies
To compare the root mean square errors (RMSE) of the proposed approach with MUSIC and ESPRIT, simulation experiments were conducted in MATLAB (Math-Works.Inc. Natick, Massachusetts, USA). In this paper, the RMSE [54] criterion was defined as where P denotes the number of signal sources, e N denotes the number of independent experiments, θ p is the pth DOA of the incident signal source, and ( ) θˆp n is an estimate of θ p in the nth Monte Carlo trial.
In the first experiment, the algorithm parameters were set as follows: there were 10 elements in the ULA, and the spacing between adjacent elements was λ = 2 d . We assumed that three incoherent signals coming from  5 ,  10 , and 15  impinged on the ULA. The SNRs of the three signals were set to be −15 dB, and the snapshot varied from 10 to 50 with 10 intervals; for the MUSIC algorithm, we set the search step as 0.1  . A total of 100 experiments were completed. Figure 3a shows the performance of the RMSE versus snapshots (signals with the same SNR). It was apparent that the proposed method outperformed the other methods. The RMSE of the proposed method became smaller as the number of snapshots increased. The principles of MUSIC and ESPRIT determine that their performance is highly dependent on the covariance matrix and the signal/noise subspace. Under the condition of a low SNR and small number of snapshots, their performance deteriorates significantly or even fails due to the poor accuracy of the covariance matrix of the array output and signal/noise subspaces. In addition, in Figure 3b, we also plot the cost function and the RMSE versus the number of iterations. The number of snapshots was fixed at 30. It was noticeable that the cost function and the RMSE decreased as the iterative process proceeded, and the method could converge within 45 iterations.  In the second experiment, the three incoherent signals were set to be 0 dB, −15 dB, and −10 dB, respectively. Other parameter settings were the same as in the first experiment. As shown in Figure 4a, the RMSE of the proposed method was smaller than that of other algorithms and was not sensitive to the change in the snapshot number. As examples of traditional characteristic subspace methods, MUSIC and ESPRIT require a large number of snapshots to construct a sampling covariance matrix, which is used to replace the data covariance matrix in subsequent signal processing. However, a small number of snapshots and a low SNR can cause a remarkable difference between the estimated covariance matrix and the real covariance matrix, which causes the DOA detection performance to degrade significantly, as shown in Figure 4.
In the second experiment, the three incoherent signals were set to be 0 dB, −15 dB, and −10 dB, respectively. Other parameter settings were the same as in the first experiment. As shown in Figure 4a, the RMSE of the proposed method was smaller than that of other algorithms and was not sensitive to the change in the snapshot number. As examples of traditional characteristic subspace methods, MUSIC and ESPRIT require a large number of snapshots to construct a sampling covariance matrix, which is used to replace the data covariance matrix in subsequent signal processing. However, a small number of snapshots and a low SNR can cause a remarkable difference between the estimated covariance matrix and the real covariance matrix, which causes the DOA detection performance to degrade significantly, as shown in Figure 4.
The proposed method transforms a multi-objective optimization problem into a single-objective optimization problem, making it possible to estimate the DOAs of multiple signal sources in a single spectral peak. The performance of the DOA detection of different SNR signals by the proposed method was much better than that of subspace-class algorithms, in a scenario with a low SNR and small snapshot number, as verified by the above simulation. The cost function and the RMSE versus the number of iterations (the number of snapshots was fixed at 30) are displayed in Figure 4b, and the method can finish implementation in 50 iterations in a situation with signals with different SNRs.

Conclusions
An augmentation of the DOA detection approach, with reconstruction of the signal subspace, is presented in this paper. During the design process, we constructed a novel cost function for the DOA through the orthogonality between the two subspaces-more precisely, signal and noise subspaces. This cost function transforms the multi-objective optimization problem in traditional subspace-class algorithms into a single-objective optimization problem, which can estimate the DOAs of multiple signal sources from a single spectral peak and avoid false alarms or missed detection in the event of a low SNR and small number of snapshots. The experimental results demonstrate that the proposed method has better DOA detection performance with a low SNR and small number of snapshots. To the best of our knowledge, this research scheme is being introduced for the first time. The proposed method transforms a multi-objective optimization problem into a singleobjective optimization problem, making it possible to estimate the DOAs of multiple signal sources in a single spectral peak. The performance of the DOA detection of different SNR signals by the proposed method was much better than that of subspace-class algorithms, in a scenario with a low SNR and small snapshot number, as verified by the above simulation. The cost function and the RMSE versus the number of iterations (the number of snapshots was fixed at 30) are displayed in Figure 4b, and the method can finish implementation in 50 iterations in a situation with signals with different SNRs.

Conclusions
An augmentation of the DOA detection approach, with reconstruction of the signal subspace, is presented in this paper. During the design process, we constructed a novel cost function for the DOA through the orthogonality between the two subspaces-more precisely, signal and noise subspaces. This cost function transforms the multi-objective optimization problem in traditional subspace-class algorithms into a single-objective optimization problem, which can estimate the DOAs of multiple signal sources from a single spectral peak and avoid false alarms or missed detection in the event of a low SNR and small number of snapshots. The experimental results demonstrate that the proposed method has better DOA detection performance with a low SNR and small number of snapshots. To the best of our knowledge, this research scheme is being introduced for the first time.
The proposed scheme puts forward a novel method for improving the performance of DOA detection and also poses a more general problem concerning the reduction of computational complexity, since the proposed method involves the QPSO optimization, which could create some potential computing overhead. Future work should consider further reducing RMSE and its applicability to more complex array configurations with the proposed method. In addition, hardware design and algorithm implementation would also be an interesting topic for discussion.