Passive Localization of Mixed Far-Field and Near-Field Sources without Estimating the Number of Sources

This paper presents a novel algorithm for the localization of mixed far-field sources (FFSs) and near-field sources (NFSs) without estimating the source number. Firstly, the algorithm decouples the direction-of-arrival (DOA) estimation from the range estimation by exploiting fourth-order spatial-temporal cumulants of the observed data. Based on the joint diagonalization structure of multiple spatial-temporal cumulant matrices, a new one-dimensional (1-D) spatial spectrum function is derived to generate the DOA estimates of both FFSs and NFSs. Then, the FFSs and NFSs are identified and the range parameters of NFSs are determined via beamforming technique. Compared with traditional mixed sources localization algorithms, the proposed algorithm avoids the performance deterioration induced by erroneous source number estimation. Furthermore, it has a higher resolution capability and improves the estimation accuracy. Computer simulations are implemented to verify the effectiveness of the proposed algorithm.


Introduction
Source localization has received considerable attention in sensor array signal processing over the past decades [1]. Most of the existing algorithms concentrate on far-field sources (FFSs), whose wavefronts are plane waves. Many high resolution algorithms have been proposed for the direction-of-arrival (DOA) estimation under the far-field assumption, such as the multiple signal classification (MUSIC) method [2], the estimation of signal parameters via rotational invariance technique (ESPRIT) [3], and their derivatives [4][5][6]. However, in many practical applications, the radiating sources may lie in the Fresnel region of the array [7], which is defined as [0.62(D 2 /λ) ½ , 2D 2 /λ] where is the wavelength of the sources and D symbolizes the array aperture. In this region, the spherical wavefronts are characterized by both DOA and range parameters [8]. Consequently, traditional FFS DOA estimation algorithms are no longer applicable for near-field sources (NFSs) localization. Fortunately, many advanced algorithms have been proposed for NFSs localization, including the 2-D MUSIC algorithm [8], the covariance approximation (CA) method [9,10], the weighted linear prediction method [11], and the rank-reduction (RARE) type algorithms [12][13][14].
Moreover, both FFSs and NFSs may coexist in many situations of interest, such as seismic exploration, electronic surveillance and speaker localization using microphone arrays. Most of the algorithms which deal with pure FFSs or pure NFSs, may fail in the scenarios of mixed sources. Accordingly, there has been an increasing interest in mixed sources localization. A two-stage MUSIC (TSMUSIC) algorithm [15] is firstly developed to localize mixed FFSs and NFSs. Based on fourth order cumulants, the TSMUSIC algorithm can successfully estimate the parameters of mixed sources. However, its computational complexity is high due to the construction of high order cumulant matrices. To relieve the computational burden, an efficient oblique projection MUSIC (OPMUSIC) algorithm is proposed in [16], which only utilizes second-order statistics. However, this algorithm suffers from severe array aperture loss. In [17], a low-complexity ESPRIT algorithm is advanced to locate the mixed far-field and near-field cyclostationary sources. Yet, array aperture loss and reduced range estimation accuracy are two slight limitations of this method. In [18], the authors extend the array aperture by utilizing a special nested sparse linear array, which can improve the estimation accuracy. Unfortunately, the range estimation suffers from spurious peaks problem in this scheme. Based on the generalized ESPRIT (GESPRIT) algorithm [19], many algorithms have been proposed to localize mixed sources [20][21][22][23]. In [22], the related far-field components are eliminated from the signal subspace to improve the estimation accuracy. But its far-field component elimination technique would bring extra estimation errors. A covariance differencing method is proposed in [20,21] to generate more reasonable classification of the source types. Based on ESPRIT-Like and polynomial rooting, an efficient mixed sources localization algorithm is presented in [23]. However, for these GESPRIT-based algorithms, the number of NFSs they can resolve is less than half of the number of array sensors [24]. In [25][26][27], the authors utilize the sparse signal recovery technique for mixed sources localization, which provide improved estimation accuracy. The algorithms in [26,27] are based on the construction of fourth order cumulant matrices and vectors, whereas the anti-diagonal elements of the second-order array covariance matrix is exploited in [25]. However, these sparse-recovery-based algorithms call for an enormous amount of computations. Besides, the regularization parameter that balances the tradeoff between Frobenius norm and norm is difficult to determine. the proposed algorithm. In Section 6, simulations are conducted to validate the performance of our method. We conclude this paper in Section 7. Throughout the paper, the complex conjugate, transpose, Hermitian transpose, pseudo inverse are denoted by (•) * , (•) T , (•) H and (•) † , respectively. Im represents a m m × identity matrix, and 0m,n is a m n × zero matrix.

Problem Formulation
Consider K independent narrowband sources (near-field or far-field) sensed by a symmetric uniform linear sensor array (ULSA), which is composed of 2M + 1 omnidirectional sensors, as shown in Figure 1. All of these sensors are placed along the y-axis, with the inter-element spacing being d. Moreover, we assume that the impinging signals are in the y-z plane. As shown in Figure 1, the azimuth angles of impinging sources become 90° in this scenario, therefore, only the elevation angles and ranges need to be estimated. However, the proposed algorithm can be easily extended to address azimuth-elevation-range (3-D) localization problem by utilizing planar arrays or placing additional sensors along the x-axis. The present signal model involves K1 sources in the far-field and the rest K2 sources in the near-field, With the array center being the phase reference point, the data observed by the mth sensor at time index t has the following form: where ( ) k s t is the kth source signal, ( ) m n t represents the mth sensor noise, mk τ is the phase shift associated with the kth source propagation time delay between the reference sensor and the mth sensor, which is of the form: where λ is the source wavelength. It is obvious that the phase shift mk τ is a nonlinear function of the Note that, for the FFSs, the second term k β is approximated by zero and the associated range parameter k r is assumed to be ∞.
In a matrix form, the array output vector x(t) can be modeled as: where s(t) is the 1 K × signal vector of the K sources: n(t) is the (2 1) 1 M + × sensor noise vector: Given the observed array data, a novel algorithm is proposed in Section 3 to localize and distinguish the mixed sources successfully, under the following hypotheses.
(1) The DOA parameters , 1, , k k K θ =  are distinct; (2) The incoming signals are mutually independent, narrowband stationary, and non-Gaussian, having nonzero fourth-order cumulants; (3) The noise is zero-mean, additive (white or color) Gaussian, and statistically independent from all impinging sources; (4) In order to avoid the phase ambiguity, the inter-sensor spacing d should be within a quarter wavelength.

Construction of the Spatial-Temporal Fourth-Order Cumulant Matrices
In this subsection, we apply the fourth-order cumulants to construct multiple spatial-temporal cumulant matrices. Due to the high degrees of freedom (DOF) available from cumulants, the resultant cumulant matrices are able to decouple the DOA estimation from the range estimation. Moreover, the array aperture size is fully utilized in these cumulant matrices so that the aperture loss problem is circumvented.
According to the definition in [34], at time lag τ , the fourth order spatial-temporal cumulant of the array outputs ( ) x t can be written as: Assuming that the source cumulants 4, ( ) sk c τ are nonzero for L different time lags l τ (1 l L ≤ ≤ ), we can formulate L fourth-order spatial-temporal cumulant matrices, which can be expressed in a matrix form:

DOA Estimation without Knowledge of Source Number
From Equation (12), it is evident that each of the L spatial-temporal cumulant matrices spans the same column space with that of the virtual steering matrix A . Therefore, with the joint diagonalization structure of the multiple cumulant matrices, we can utilize these matrices simultaneously to identify the column space of A and generate the DOA estimates of the mixed sources.
For the kth source, one can always find a vector which is orthogonal to the range space spanned by the virtual steering vectors except ( ) (13) and the following notation can be equally used: Substituting Equation (14) into Equation (12) leads to: where θ is the DOA parameter of interest, C are the nuisance parameters. In order to avoid the trivial solution {b=0, d=0}, the constraint 1 = d is added. From Equation (16), it is obvious that the computational complexity of this estimator is prohibitive when b and d are unknown. Therefore, to decrease the computational load, we first decouple the DOA parameter θ from other parameters, where only 1-D spectral search is required.
The objective function in Equation (16) can be expanded as: Let: In view of that  (17) can be rewritten as: In order to decouple the nuisance parameter b from the objective function, we set the partial derivative of ( , , ) which yields: Substituting Equation (22) into Equation (20), the optimization problem without the nuisance parameter b is of the following form: The last equation holds when 1 = d v is the eigenvector corresponding to the maximum eigenvalue of ( ) , which is 1 λ . Therefore, the estimation of DOA is further simplified as: Herein, 'max eig' stands for the maximum eigenvalue of a square matrix. Consequently, the DOA estimates of the mixed sources can be obtained from the highest peaks of the following 1-D spectral function:

Range Estimation and Source Classification
In order to circumvent the source number enumeration, the minimum variance distortionless response (MVDR) beamformer [35,36] is utilized in this subsection for the range estimation problem. The output power spectrum of the MVDR beamformer is: where x x is the second-order array covariance matrix. With the DOA estimates obtained in the previous subsection, the 2-D spectrum search over the DOA-range plane could be reduced to 1-D search on the range domain. Therefore, by substituting the estimated DOA ˆk θ back into Equation (28), the range k r can be obtained as: Since the NFSs lie in the Fresnel region, r is searched over [0.62(D 2 /λ) ½ , 2D 2 /λ]. If the estimation ˆk r is in the Fresnel region, the corresponding source is identified as a near-field one; otherwise, it is determined to be in the far-field. Note that, in this procedure, the DOA estimation and range estimation are automatically paired without any additional operation.

Implementation of the Proposed Algorithm
Based on the above analyses, the implementation of the proposed algorithm is summarized as follows: Step 1. Calculate the L fourth-order spatial-temporal matrices Step 2. Construct the matrices C and ( ) θ Q from Equations (18) and (19); Step 3. Use Equation (27) to plot the DOA spectrum function ( ) P θ , and find the DOA estimates; Step 4. Find the range estimates using Equation (29), and classify the types of the sources;

Analogy between the Proposed Algorithm and TSMUSIC
In this subsection, we analyze the relationship between the proposed algorithm and TSMUSIC in the framework of subspace decomposition. If we take only one single spatial-temporal cumulant matrix  (18) and (19) will be simplified as: and consequently, can be reformulated as: Accordingly, the spectrum function in Equation (27) becomes: The eigenvalue decomposition (EVD) of ( ) where s Λ is the diagonal matrix containing K large eigenvalues. Es is the ( ) which is identical to Equation (35). Therefore, we come to an interesting conclusion that the DOA estimator for TSMUSIC is a special case of our method, if only one spatial-temporal cumulant matrix is utilized. However, our method does not require the estimation of the source number, which makes it practically attractive.

Discussion
In this section, we compare the proposed algorithm with two newly developed mixed sources localization algorithms: TSMUSIC [15] and OPMUSIC [16]. We discuss all of these three methods from the following aspects: (1) Computational complexity: Concerning the computational complexity, we only compare the major multiplications involved in the construction of the statistical matrices, eigenvalue decomposition (EVD) and spectral search. The searching steps for the DOA parameter and the range parameter are denoted as θ Δ and r Δ . Let T and 2M + 1 symbolize the snapshot number and sensor number, respectively. TSMUSIC requires computing two fourth-order matrices with dimension ( ) ( )  . Therefore, the computational complexity of the proposed algorithm is higher than that of TSMUSIC and OPMUSIC. However, it is important to note that the proposed algorithm does not require the knowledge of the source number, which is highly desirable for practical applications where detection of the source number is generally a very difficult task.   (14). In other words, the maximum number of resolvable sources for the proposed algorithm is 2M + 1 with the same sensor array configuration. As a result, our method outperforms TSMUSIC and OPMUSIC in resolving more sources.

Simulation Results
In this section, several numerical simulations are conducted to validate the performance of the proposed algorithm relative to TSMUSIC [15] and OPMUSIC [16]. In the following simulations, we consider an symmetric ULSA composed of 2M + 1 = 9 (M = 4) elements with . The impinging signals are narrowband, noncoherent, equi-power and with the non-Gaussian form . For the proposed algorithm, L = 5 spatial-temporal cumulant matrices at the first five time lags are considered. Moreover, it is always assumed that the source number is correctly estimated for TSMUSIC and OPMUSIC. The performance is measured in terms of spatial spectrum, probability of resolution (PR) and root mean-square error (RMSE). The PR is defined as the ratio between the number of successful resolution and the total number of Monte Carlo runs. If the DOA estimates of two sources , they are said to be successfully resolved. Herein, signifies the angular separation between the two sources. The RMSE is defined as: where stands for the parameters such as and , denotes the estimation of in the nth trial, and is the number of Monte Carlo runs.
In the first experiment, we investigate the maximum number of resolvable sources for the proposed algorithm. We consider the mixed sources scenario that four FFSs are located at , , , , and five NFSs are placed at , , , , , respectively. The snapshot number and the SNR equal to 200 and 20 dB, respectively. Figure 2 indicates the DOA spectra from 10 independent realizations. It is evident that, with the 9-element ULSA, there are nine clearly discernible peaks corresponding to these mixed sources, which is in agreement with the analysis in Section 5. number and the SNR equal to 500 and 5 dB, respectively. Ten independent trials of DOA spectrum estimates using the three methods are realized as is shown in Figure 3. It is obvious that, for the proposed method, there are two distinct peaks corresponding to the actual DOAs of the two mixed sources. However, TSMUSIC and OPMUSIC fails to distinguish the two closely spaced angles 1 θ and 2 θ . This is because that, for closely spaced sources, it is very difficult for the subspace-based algorithms to correctly identify the signal subspace and noise subspace with low SNR, while the proposed method does not need any decomposition and identification of the subspaces.
In the third experiment, the PR of the three algorithms versus SNR is explored. Consider two uncorrelated equi-power sources locating at SNR, 500 independent Monte Carlo trials are performed. Figure 4 illustrates the PR of the three algorithms as a function of SNR. It is obvious that the proposed method has the best resolution ability since it reaches a full PR at the lowest SNR threshold (1dB). For TSMUSIC and OPMUSIC, the SNR thresholds of full PR are 10 dB and 25 dB, respectively, which are much higher than that of our method. This is due to the fact that TSMUSIC and OPMUSIC suffer from the leakage between the signal subspace and noise subspace when SNR is low, while the subspace decomposition is not required in our algorithm, which promotes the resolution ability at low SNRs. Moreover, OPMUSIC has the worst resolution ability since its DOA estimator has a half aperture loss, while the array size is fully utilized in TSMUSIC and our method.   In the fourth experiment, we study the estimation accuracy of the three algorithms as a function of SNR. A more general mixed sources scenario is considered, where two FFSs are located at  Figures 5-7 illustrate the RMSEs of DOA and range estimates using the three algorithms. It can be seen that the estimation accuracy (both DOA and range) of the proposed algorithm outperforms that of TSMUSIC and OPMUSIC, especially in the low SNR region. This can be also explained from the fact that incorrect identifications of signal subspaces are likely occur at low SNRs, which would deteriorate the performance of TSMUSIC and OPMUSIC. At the high SNR region, the DOA estimation accuracy of TSMUSIC approaches to that of our method, which is in agreement with the analysis presented in Section 4. OPMUSIC has the worst DOA estimation accuracy because of its half aperture loss. As for the range estimation, in the low SNR region, the proposed algorithm is still superior to TSMUSIC and OPMUSIC, which is mainly a result of the propagation error from the previous DOA estimation stage. In the high SNR region, however, OPMUSIC and our method have almost the same performance. Additionally, TSMUSIC has the worst range estimation performance since it only utilizes the information of a single eigenvector with the smallest eigenvalue to obtain the range estimates. As is described in [15], the DOA estimates and the range estimates are obtained from two different stages in TSMUSIC. In the first stage, a special cumulant matrix is constructed and decomposed. DOAs can be obtained via high resolution MUSIC algorithm. In the second stage, another particular cumulant matrix is constructed to avoid the estimation failure problem. However, it only utilizes the information of a single eigenvector with the smallest eigenvalue to obtain the range estimates in this stage. Therefore, TSMUSIC has a good performance in bearing estimation but a worst performance in range estimation in the simulations.    In the fifth experiment, the same parameters as the third experiment are adopted except that the SNR is set equal to 15 dB, and the number of snapshots varies from 10 to 1000. At each snapshot number, 500 independent Monte Carlo trials are performed. The RMSEs of DOA and range estimation are shown in Figure 8-10. It is obvious that the RMSEs (both DOA and range) of the three algorithms decrease monotonically as the snapshot number increases. This is due to the fact that a larger sampling number will produce better estimate of the cumulant matrices and the covariance matrices for stationary data. Furthermore, when the snapshot number is large enough, the proposed algorithm slightly outperforms TSMUSIC and OPMUSIC. This is because that: (1) the array aperture is fully utilized in this algorithm; (2) both the temporal and spatial structures of the observed data are utilized in our method, which also contribute to the improvement of estimation accuracy. However, when the sample size is small (i.e., snapshot number = 10), the DOA estimation accuracy of our method is a little bit inferior to that of TSMUSIC and OPMUSIC. This can be explained that the estimation of the spatial-temporal cumulant matrices is not reliable when the sample size is insufficient.    In the last experiment, we compare the computational complexity of the three algorithms. Suppose that there are 2M + 1 = 9 sensors, and that the searching steps θ Δ and r Δ are set as 0.1° and 0.05λ , respectively. The number of snapshots varies from 10 to 1000. Furthermore, we define K = 2 and L = 5. Figure 11 illustrates the computational burden of these algorithms as a function of the snapshot number. It is clear that the proposed algorithm has slightly higher complexity than TSMUSIC, while OPMUSIC is most efficient in computational complexity since it only involves second-order statistics. Figure 11. Computational complexity of the three methods versus snapshot number.

Conclusions
In this paper, a novel algorithm has been proposed for the localization of mixed FFSs and NFSs. Based on the joint diagonalization structure of multiple fourth-order spatial-temporal cumulant matrices and the beamforming technique, the DOA and range parameters of FFSs and NFSs have been determined respectively. Compared with some existing mixed source localization methods, the new method has the main advantage that it does not require the information of the source number. Such an advantage is practically attractive since mixed source enumeration is typically a very difficult problem. In addition, we have found that our method is a generalization of the TSMUSIC algorithm in [15], which only utilize one single spatial-temporal cumulant matrix. Moreover, the proposed algorithm avoids multidimensional search, aperture loss and parameter matching. According to the simulation results, the proposed method outperforms TSMUSIC and OPMUSIC in both resolution ability and estimation accuracy, especially when the SNR is low.
It is worth mentioning that there are two slight limitations of the proposed algorithm. First, the construction of cumulant matrices leads to high computational load. Second, it is limited to localize mixed sources only in the 2-D domain (elevation angle and range). To put the proposed algorithm into further applications, several steps of the future work can be shown as follows: 1. To reduce the computational complexity, second-order statistics and real-valued transformation may be applied. 2. To extend the proposed algorithm for joint azimuth-elevation-range (3-D) estimation problem, we may employ planar arrays, such as cross array, rectangular array and circular array.