Mixed Far-Field and Near-Field Source Localization Algorithm via Sparse Subarrays

Based on a dual-size shift invariance sparse linear array, this paper presents a novel algorithm for the localization of mixed far-field and near-field sources. First, by constructing a cumulant matrix with only direction-of-arrival (DOA) information, the proposed algorithm decouples the DOA estimation from the range estimation. The cumulant-domain quarter-wavelength invariance yields unambiguous estimates of DOAs, which are then used as coarse references to disambiguate the phase ambiguities in fine estimates induced from the larger spatial invariance. Then, based on the estimated DOAs, another cumulant matrix is derived and decoupled to generate unambiguous and cyclically ambiguous estimates of range parameter. According to the coarse range estimation, the types of sources can be identified and the unambiguous fine range estimates of NF sources are obtained after disambiguation. Compared with some existing algorithms, the proposed algorithm enjoys extended array aperture and higher estimation accuracy. Simulation results are given to validate the performance of the proposed algorithm.


Introduction
In recent years, passive source localization has become a key topic in array signal processing [1].Various localization algorithms have been proposed for far-field (FF) source, whose wavefronts are plane waves, such as the multiple signal classification (MUSIC) method [2], the estimation of signal parameters via rotational invariance technique (ESPRIT) [3], and their derivatives.Nevertheless, when the radiating sources are located in near-field (NF) source, whose wavefronts are spherical waves, both the DOA and the range parameters should be determined to localize these radiating sources.As a result, traditional FF DOA estimation algorithms are no longer applicable for NF source localization.Fortunately, many advanced methods have been presented under the NF assumption, including the 2-D MUSIC algorithm [4], the high-order ESPRIT algorithm [5,6], the covariance approximation (CA) method [7,8], the weighted linear prediction method [9], the generalized ESPRIT algorithm [10,11], and the signal reconstruction for nearfield source localization [12].In addition, we also proposed a novel method of passive localization for near-field noncircular sources [13].
However, both FF and NF sources may coexist in many interested situations such as speaker localization using microphone arrays, seismic exploration, and electronic surveillance.Most of the algorithms, which deal with pure NF or pure FF sources, may fail in the scenarios of mixed sources.
Recently, mixed source localization problem has been an important research topic in array signal processing [14][15][16][17][18][19].A two-stage MUSIC (TSMUSIC) algorithm [14] was first advanced to localize mixed FF and NF sources.Based on fourth-order cumulant, 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 and the spectral search.To relieve the computational burden, an efficient MUSIC-based algorithm is proposed in [15], which only utilizes secondorder statistics.Unfortunately, this algorithm suffers from severe array aperture loss and 1-D spectral search in both DOA and range estimation.Based on [10], Liu and Sun presented another GESPRIT-based algorithm to alleviate the array aperture and obtain a reasonable classification result [16].Our early works on mixed localization focused on unknown source numbers [17] and unknown mutual coupling [18].
As is well known, the estimation accuracy is directly correlated with the array aperture size that a larger array would produce more precise estimates.However, most of the existing methods limit the array element spacing to be within a quarter wavelength to avoid DOA ambiguity.Recently, a mixed-order MUSIC (MOMUSIC) algorithm [19] was proposed to extend the array aperture by a special nested sparse linear array (SLA), which can improve the estimation accuracy.
In this paper, a novel fourth-order cumulant-based dualsize shift invariance (CDSSI) algorithm is presented to solve the mixed source localization problem.Our technique utilizes a SLA of the dual-size spatial invariance method [20] which was designed by M. D. Zoltowski and K. T. Wong.Furthermore, they also proposed MUSIC/MODE null spectrum disambiguation algorithm [21] aiming to realize more accurate disambiguation.Unlike most of the existing algorithms, our technique utilizes a SLA of dual-size spatial invariance.As a result, the novel method enjoys significant promotion in DOA and range estimation accuracy by extending the intersubarray spacing.Furthermore, the proposed method avoids any 1-D or 2-D spectral searching and therefore has lower computational complexity.
The rest of the paper is organized as follows.In Section 2, the mixed FF and NF signal model based on SLA is presented.The proposed CDSSI algorithm is described in Section 3. In Section 4, we compare CDSSI with some recently developed ones, and computer simulations are conducted in Section 5 to validate the performance of the proposed algorithms.Finally, we conclude this paper in Section 6.

Data Model
Suppose that K-independent narrowband sources (FF and NF) impinge upon a symmetric SLA, as shown in Figure 1.The array has a total number of 2M + 1 sensors, which are composed of 3 subarrays and each subarray contains 2M u + 1 array sensors, with M and M u being positive integers.This SLA configuration can be considered as a particularization of the sparse rectangular dual-size spatial invariance array in [20].
In Figure 1, the intersensor spacing is d, and the intersubarray spacing is d s = d ⋅ Δ, where d s ≫ d.The sensor position vector P would have the following form, if we take d as the length unit.

Source
Figure 1: Sparse linear array geometry.

International Journal of Antennas and Propagation
Let the array center be the phase reference point; the output of the mth sensor can be approximated as [5] where θ k represents the DOA of the kth source, r k stands for the distance between the kth source and the reference sensor, s k t symbolizes the kth source signal, n m t denotes the additive Gaussian noise, and L is the number of snapshots.It should be noted that, for the NF source scenario, the range parameter lies in the Fresnel region 0 62 D 3 /λ, 2D 2 /λ [22], with D symbolizing the array aperture.For the FF source scenario, the range parameter approaches to ∞ and the associated parameter ϕ k becomes 0 [14].In a matrix form, the array data can be written as where In the above equations, A represents the array steering matrix, a ω k , ϕ k denotes the array steering vector for the kth source, S t symbolizes the source signal vector, and N t is the noise vector.
Given the array data X t , a novel algorithm is proposed in Section 3 to obtain high performance in localizing and distinguishing the mixed sources successfully, under the following hypotheses.
(i) The incoming signals are mutually independent, narrowband stationary, and non-Gaussian, as well as nonzero kurtosis.
(ii) The DOAs of all source signals differ from each other.
(iii) The noise is zero-mean, additive (white or color) Gaussian, and statistically independent from all impinging sources.
(iv) In order to avoid the phase ambiguity, the intersensor spacing d should be within a quarter wavelength.

Algorithm Development
As is shown in Figure 2, the proposed algorithm has two main stages.In the first stage, by constructing a fourth-order cumulant matrix C 1 containing only the DOA information, we can decouple the DOA estimation from the range estimation.After the eigendecomposition of C 1 , dual-size shift invariance ESPRIT [20] could be applied to generate coarse (unambiguous) and fine (ambiguous) DOA estimates of both NF and FF sources.In the second stage, another fourth-order cumulant matrix C 2 needs to be constructed to overcome the rank-deficient phenomenon described in [14].Moreover, the virtual steering matrix of C 2 can be decoupled into two parts and both the coarse and fine range estimates of these sources could be obtained.In both of the two stages, disambiguation is important for the final accurate and unambiguous localization of mixed FF and NF sources.The proposed CDSSI algorithm is described in detail in the following subsections.

DOA Estimation
3.1.1.Cumulant Matrix for the DOA Estimation.Similar to the TSMUSIC algorithm in [14], we first construct a fourth-order cumulant matrix to estimate the DOAs of the radiating sources in this subsection.However, there are some differences between the two algorithms.For example, (1) the uniform linear array (ULA) configuration in TSMUSIC is generalized to a SLA in Figure 1, which promotes the estimation accuracy effectively and (2) the proposed method utilizes ESPRIT to generate the DOA estimation, which avoids the spectral search procedure in TSMUSIC and therefore reduces the computational complexity.
According to the definition in [23], the fourth-order cumulant of the array outputs x m t , x n t , x i t , and x j t can be written as where symbolizes the kurtosis of the kth source.
We can construct a 2M + 1 × 2M + 1 cumulant matrix C 1 , with its m, n th element being Note that C 1 can be written in a matrix form where From (1), B can be written as where International Journal of Antennas and Propagation herein, A sub is the virtual steering matrix for each subarray in the cumulant domain, with the following form: Therefore, C 1 can be expressed as The reason to rewrite it in this manner is that it has a similar form of dual-size invariance and therefore, ESPRIT could be applied to yield coarse and fine estimates of DOAs.
3.1.2.Fine DOA Estimates with Ambiguities.We firstly estimate the DOAs of the mixed FF and NF sources.By taking an eigendecomposition of (17), we have where Λ S and Λ N are diagonal matrices containing K large eigenvalues and 2M + 1 − K small eigenvalues, respectively.U N is the 2M + 1 × 2M + 1 − K eigenvector matrix spanning the noise subspace of C 1 .U S is the 2M + 1 × K signal subspace eigenvector matrix of C 1 , which can be expressed as where T is a nonsingular matrix.
To obtain high-accuracy DOA estimates, we may form two matrices that are related by the extended intersubarray spacing d s along the x-axis.As is illustrated in Figure 3, by taking the first and the last 2 2M u + 1 rows of U S , we have the following two matrices.

20
where D 1 and D 2 are the first and the last 2 rows of D s and Φ f = diag e j2Δω 1 , … , e j2Δω K .Since U S1 is contributed by the cumulant domain data from all 2 2M u + 1 sensors at subar ray −1 and subarray 0 and U S2 is contributed by the cumulant domain data from all 2 2M u + 1 sensors at subar ray 0 and subarray 1 , fine DOA estimates can be found from Φ f due to the extended invariance relationship.
From (20), we have e j4πd s sin θ k /λ , correspond to the eigenvalues of Ψ f .As d s ≫ λ/4 and sin θ k ≤ 1, we can find a series of ambiguous DOA estimates where l θ is an integer and ⋅ and ⋅ represent the ceil (round toward positive infinity) and floor (round toward negative infinity) operation in MATLAB, respectively.∠ ⋅ returns the phase angle of the operand, which lies between −π, π .
3.1.3.Coarse DOA Estimates without Ambiguities.For the purpose of disambiguation, unambiguous but coarse DOA estimates must be obtained as references to the ambiguous fine estimates in the previous subsection.From Figure 4, by extracting the cumulant domain information of the first and the last 2M u sensors in each subarray, the quarterwavelength spatial invariance of the array geometry can be exploited to generate unambiguous coarse DOA estimates.This can be easily done by defining a permutation matrix where e i signifies the ith column of the identity matrix I 2M u +1 .From ( 19) and ( 24), we have where A 1 and A 2 are the first and the last 2M u rows of A sub , respectively.Φ c = diag e j2ω 1 , … , e j2ω K .Since U P1 is contributed by the cumulant domain data from the first 2M u sensors in each subarray and U P2 is contributed by the cumulant domain data from the last 2M u sensors in each subarray, unambiguous DOA estimates with quarter-wavelength spatial invariance can be found from Φ c .From (26), we have U P2 = U P1 T −1 Φ c T, which means that there is a rotational invariance between U P1 and U P2 , that is, Therefore, the diagonal elements of Φ c , Φ c k,k = e j2ω k = e j4πd sin θ k /λ , correspond to the eigenvalues of Ψ c .The unambiguous coarse DOA estimates of K-mixed FF and NF sources are given by θref k and θk l θ may be paired by the method in [20].3.1.4.Disambiguation.The coarse estimates θref k could be used as references to disambiguate the cyclic ambiguities in θk l θ .Mathematically, we obtain the disambiguated angle estimates as where the estimate of l 0 θ is given by Note that l 0 θ is solved through searching over a few discrete points and the searching range is given by ( 22).

Range Estimation
3.2.1.Cumulant Matrix for the Range Estimation.In order to derive the range estimates, another fourth-order cumulant matrix needs to be constructed to overcome the rankdeficient phenomenon described in [14].The virtual steering matrix for range estimation in [5] has the following form If the kth source is in FF, the electrical angle ϕ k will become 0 and the kth column of A will be 1, … , 1 T .When the number of FF sources is more than one, A will drop rank and thus the algorithm will fail to localize radiating sources.
In this section, a new fourth-order cumulant matrix C 2 is defined, with its m, n th element being Note that C 2 can be expressed as where A is the steering matrix in (6) and C 4s is the kurtosis diagonal matrix of sources as defined in (12).
When the kth source is in FF, the corresponding column vector a ω k , ϕ k becomes a ω k , ϕ k = e jω k p −M , … , e jω k p 0 , … , e jω k p M T 33 Therefore, A will still be of full column rank when there are multiple FF sources.6 International Journal of Antennas and Propagation

Unambiguous and Ambiguous Range Estimations.
From ( 1) and ( 7), a ω k , ϕ k can be decoupled as where V ω k is a 2M + 1 × M + 1 matrix which only contains the DOA information and g ϕ k is a M + 1 × 1 column vector which only depends on ϕ k .By taking an eigendecomposition of C 2 , we can obtain a 2M + 1 × K signal subspace matrix E S and a 2M + 1 × 2M + 1 − K noise subspace matrix E N , where E S is composed of the K principle eigenvectors of C 2 and E N contains the remaining 2M + 1 − K eigenvectors.
According to [4], the 2-D MUSIC spectrum for DOA and range estimation is given by Substituting the disambiguated DOA estimates θk k = 1, 2, … , K into the above equation, the 2-D spectrum search can be reduced to the K 1-D ones.Therefore, the estimates of Actually, (37) implies that g ϕ k is the eigenvector associated with the minimum eigenvalue of the Hermitian matrix V H ωk E N E H N V ωk .Through the eigendecomposition of V H ωk E N E H N V ωk , we can obtain an estimation of g ϕ k .
In order to generate the coarse range estimates unambiguously, ϕ k can be estimated as where g ϕ k n represents the nth element in g ϕ k .From (3), the coarse range estimate is given by where θk is the disambiguated DOA estimation.According to the definition of the Fresnel region, NF sources are located in the range of 0 62 D 3 /λ , 2D 2 /λ .As a result, when rref k is in this region, the corresponding source is classified as a NF one; otherwise, it is regarded as a FF one.
Similarly, through the extended aperture size between each subarray, fine range estimates with ambiguity can be generated.From (35), we have Therefore, we can find a series of ambiguous range estimates 3.2.3.Disambiguation.In an analogous manner to that of Section 3.1.4,the coarse estimates rref k could be used as references to disambiguate the cyclic ambiguities in rk l r .Mathematically, the disambiguated range estimates are where the estimate of l 0 r is given by Note that l 0 r is solved through searching over a few discrete points, and the searching range is given by (42).

Discussion
In this section, we compare CDSSI with some recently developed mixed source localization algorithms, including TSMUSIC [14], the generalized ESPRIT-based (GESPRIT) algorithm [16], and MOMUSIC [19].All of the four methods are analysed from the following aspects: 7 International Journal of Antennas and Propagation 4.1.Array Aperture.Both TSMUSIC and GESPRIT employ an ULA which requires the intersensor spacing to be within a quarter wavelength.Therefore, with the same sensor number, SLA enjoys extended array aperture size, producing better estimation accuracy for the proposed algorithm.MOMUSIC also utilizes a special nested SLA in the development of the algorithm.However, according to the array model described in [19], the array aperture of MOMUSIC equals to when M 1 = 2 and M 2 = 2 (9 sensors in all), while the array aperture of CDSSI, with the same number of sensors, could be extended by spacing subarrays apart at tens of quarter wavelengths or more [20].

Maximum Number of Sources
That Can Be Resolved.In CDSSI, the coarse and fine DOA estimates are derived from two 6M u + K matrices and two 2 2M u + 1 × K matrices, respectively; therefore, the maximum number of sources that can be resolved is 4M u + 1.For MOMUSIC algorithm, a virtual ULA with 2 M 1 + 1 M 2 + 1 + 1 sensors can be constructed.However, it has a half aperture loss in the formation of a Toeplitz matrix, and therefore, it can resolve M 1 + 1 M 2 + 1 sources at most.With a 2M + 1 element ULA, TSMUSIC is capable of resolving up to 2M sources, while GESPRIT can resolve a maximum number of 2M − 1 sources.

Computational Complexity.
Only the major computation load is considered in this comparison, including construction of the cumulant matrices, eigenvalue decomposition (EVD), and spectral search.The searching steps for the angle parameter θ ∈ −90 °, 90 °and the range parameter r ∈ 0 62 D 3 /λ , 2D 2 /λ are denoted as θ Δ and r Δ .Let L and 2M + 1 symbolize the snapshot number and sensor number, respectively.For MOMUSIC, we assume that M 1 = M 2 = M, and therefore, the sensor number is 4M + 1. TSMUSIC involves computing two fourth-order matrices with dimensions 2M + 1 × 2M + 1 and 4M + 1 × 4M + 1 , two EVDs, and spectral search for DOA estimation.The GESPRIT algorithm requires the construction of two second-order covariance matrices, two EVDs, and spectral search for DOA and range estimation.MOMUSIC requires the construction of a fourth-order cumulant matrix and a second-order covariance matrix, two EVDs, and spectral search for DOA and range estimation.CDSSI does not need any spectral search operation but requires two fourth-order cumulant matrices with dimension . The computational complexity of the four methods is listed in Table 1.
Compared with the second-order-based methods, CDSSI requires more computations in constructing the fourth-order cumulant matrix and EVDs.However, it avoids any 1-D or 2-D complicated spectral search.Since the search steps need to be dense enough for the spectral search-based algorithms to approach their theoretical bounds, the computation load of these algorithms will in turn increase dramatically.

Performance in Correlated
Noise.Both TSMUSIC and the proposed algorithm are capable of suppressing additive Gaussian colored noise since they apply fourth-order cumulant in the whole estimation procedure, while GESPRIT and MOMUSIC, which rely on the second-order statistics, will degrade in the presence of spatially correlated noise.4.5.Arc Length and First-Order Curvature.Furthermore, knowledge of the "manifold shape" not only is essential for the investigation of ambiguities and assessment of the detection-resolution capabilities of an array but it may also prove useful in developing new and more effective methods for its search process.We studied on the two array manifold properties, namely, arc length and first-order curvature, and analyse the accuracy and resolution capabilities of mixed sources [24].
It has been shown that the response of sparse array towards a far-field/near-field source emitting narrowband spherical wavefront from azimuth θ and range r can be written as , 46 where the associated parameter is ω k = − 2πd/λ sin θ, ϕ k = πd 2 /λr k cos 2 θ, and P is the sensor position vector.
For the far-field source scenario, the range parameter approaches to ∞ and the associated parameter ϕ k becomes 0.
The rate of change of arc length ds/dθ and first-order curvature k 1 of the θ parameter curves are as follows:

48
where a implies differentiation with respect to θ and the arc length s is defined as International Journal of Antennas and Propagation

International Journal of Antennas and Propagation
The Cramer-Rao lower bound under the assumption of spatially and temporally uncorrelated large number of snapshots L ≫ 1 may be written as , 50 where the unit-norm tangent vector u 1 θ, r to the array manifold has been substituted for a θ, r /s θ, r .The expression of the asymptotic (L ≫ 1) variance for a single emitter of unit power can be written as The variance asymptotically approaches the CRB for high SNR or large N and has a similar dependence on s θ, r .
Figures 5, 6, and 7 show the Cramer-Rao lower bound capabilities of our proposed sparse array.Figures 5 and 6 indicate the variation of theoretical CRB with respect to the change of r. Figure 7 shows the theoretical CRB with respect to the change of snapshot.From Figures 5 and 6, we can see that the smaller the angle with the normal direction of array, the less the estimation error and when the range between the array and source is enlarged, the theoretical estimation performance becomes worse.And we can conclude that more sample data can improve estimation performance from Figure 7.
In addition, the angle and range estimation accuracy is not only related to the intersubarray spacing of sparse but also related to the possibility of disambiguation.As the intersubarray spacing increases, the possibility of wrong disambiguation increases simultaneously.The MIE method [25] can be used to predict MSE performance; the angle and range estimation can be represented as

52
where u and û denote the estimation value and true value of angle and range parameters, P interval error is the possibility of right disambiguation, and P no interval error is the possibility of wrong disambiguation.

Simulation Results
In this section, numerical simulations are conducted to validate the performance of the proposed algorithm relative In the first experiment, we consider two equipower sources that are located at θ 1 = 10 °, r 1 = 15λ and θ 2 = −10 °, r 2 = +∞ , that is, a mixed FF and NF scenario.The number of snapshots equals to 500, and the SNR varies from −10 dB to 20 dB in steps of 3 dB.The RMSEs of the four algorithms as a function of SNR are plotted in Figure 8. From these figures, one can observe that the coarse DOA and range estimates of the proposed algorithm have higher RMSEs due to the quarter-wavelength spatial shift invariance.However, after disambiguation, the fine estimates have superior estimation accuracy than those of the other three algorithms.Moreover, RMSEs of the DOA and range estimates decrease as the SNR increases.
In the second experiment, we investigate the RMSEs of the four algorithms with the variation of the number of snapshots.The parameter settings are the same as those of the first experiment except that SNR is set equal to 10 dB and the number of snapshots varies from 100 to 10,000.From Figure 9, it is obvious that, as a result of the extended aperture, CDSSI outperforms the other three algorithms in DOA and range estimation accuracy for all snapshot numbers.In addition, both the DOA and range estimation performance of all four algorithms improves as the snapshot number increases.This is because that larger sample support will produce better estimate of the covariance matrix for stationary data.
In the third experiment, the scenario of two NF sources is investigated, with the source location parameters being θ 1 = −10 °, r 1 = 10λ and θ 1 = 10 °, r 1 = 20λ .The snapshot number is fixed at 500 and the SNR varies from −10 dB to 30 dB in steps of 5 dB. Figure 10 leads to a similar conclusion as in the first experiment that the proposed algorithm achieves the best performance owing to its extended aperture.Additionally, from the second figure, one can observe that the range estimation accuracy of the first source, which is closer to the array, is better than that of the second source.This result is consistent with the theoretical analysis developed in [6].
In the fourth experiment, we study the dependence of DOA and range estimation accuracy upon the angular gap between two NF sources.θ 1 varies form −80 °to 80 °in steps of 10 °, while θ 2 = 0 °, r 1 = 10λ, and r 2 = 20λ, with the snapshot International Journal of Antennas and Propagation From the first figure, it is obvious that the DOA estimation performance is insensitive to the change of the range parameter, since the DOA estimation is decoupled with the range estimation in all these algorithms.However, the GESPRIT algorithm behaves abnormally when the range of the first source approaches to 20λ (the range parameter of the second source).This is because when the two sources are symmetrical with respect to the broadside, that is, θ 1 = −θ 2 and r 1 = r 2 , the GESPRIT algorithm will generate image sources, which are misidentified as real ones.Moreover, the second figure reveals that (1) the range estimation performance of the second source is hardly affected by the range variation of the first source and (2) the range estimation accuracy of the first source (closer to the array) is superior to that of the second one, which corroborates the theoretical analysis in [6].Finally, from Figure 12, it is clear that the proposed algorithm outperforms the other three methods in both DOA and range estimation for all range parameters.
In the last experiment, two FF sources are considered with their location parameters being θ 1 = 10 °, r 1 = +∞ and θ 2 = 60 °, r 2 = +∞ .The snapshot number is fixed at 500 and the SNR varies from −10 dB to 20 dB in steps of 3 dB.Note that the intersubarray spacing d s is extended from 20d to 40d without introducing additional sensors.The RMSEs of the DOA estimates with the variation of SNR are shown in Figure 13.From this figure, it is seen that all of the four algorithms are still effective in the multiple FF source scenarios and CDSSI outperforms the other three in the performance of DOAs estimation.Also, one can observe that as the intersubarray spacing becomes larger, the accuracy of the proposed method turns better compared with the results demonstrated in Figure 2.

Conclusion
In this paper, an efficient and high-performance algorithm is proposed for the mixed far-field and near-field source localization problems.Based on a sparse linear array of dual-size spatial invariance, the proposed algorithm can offer enhanced accuracy due to the extended aperture size.Moreover, the proposed method has lower computational complexity because it does not require any 1-D or 2-D spectral search.According to the simulations, the proposed algorithm outperforms the conventional ones in the performance of both angle and range estimation.

Figure 2 :
Figure 2: Flow graph of the CDSSI algorithm.

3
International Journal of Antennas and Propagation

Figure 3 :
Figure 3: Construction of 2 subspace matrix elements for fine DOA estimation.

Figure 4 :
Figure 4: Construction of 2 subspace matrix elements for coarse DOA estimation.

Table 1 :
Computational complexity of different algorithms.