Sparse Reconstruction Based Robust Near-Field Source Localization Algorithm †

Non-Gaussian impulsive noise widely exists in the real world, this paper takes the α-stable distribution as the mathematical model of non-Gaussian impulsive noise and works on the joint direction-of-arrival (DOA) and range estimation problem of near-field signals in impulsive noise environment. Since the conventional algorithms based on the classical second order correlation statistics degenerate severely in the impulsive noise environment, this paper adopts two robust correlations, the fractional lower order correlation (FLOC) and the nonlinear transform correlation (NTC), and presents two related near-field localization algorithms. In our proposed algorithms, by exploring the symmetrical characteristic of the array, we construct the robust far-field approximate correlation vector in relation with the DOA only, which allows for bearing estimation based on the sparse reconstruction. With the estimated bearing, the range can consequently be obtained by the sparse reconstruction of the output of a virtual array. The proposed algorithms have the merits of good noise suppression ability, and their effectiveness is demonstrated by the computer simulation results.


Introduction
The source localization problem in array signal processing is an important problem which has a wide range of applications. For example: related topics in radar, sonar, wireless communications, and seismology is to determine the location of the radiation (or reflection) sources by passive sensor array [1][2][3][4]. According to the distance of the source from the array, the localization problem is divided into far-field source localization problem and near-field source localization problem. When the sources are far away from the array, the localization problem in this case is far-field source localization and the signals transmitted by the sources arrive at the array in the form of a plane wave. At this point, the direction of arrival (DOA) needs to be estimated. When the distance between the sources and the array is within the Fresnel region, the localization problem at this time is near-field source localization and the signals transmitted by the sources arrive at the array in the form of a spherical wave. At this time, the DOA and range must be determined. Compared with the DOA estimation in far-field, the position parameter estimation problem of near-field source is relatively complex.
In recent years, near-field source localization problem has caused the concern in the field of signal processing, and some scholars have started to research on this problem. Huang et al. [5] extended the one-dimensional MUSIC method used in the far-field to the two-dimensional MUSIC method to be used in the near-field. The biggest limitation of this method was that it required a two-dimensional spectral peak search and required a large amount of computation. Starer et al. [6] proposed a path-following approach which transformed the two-dimensional search problem into multiple one-dimensional search problems, and finally solved the problem by iterative method. Lee et al. [7] proposed a modified path-following method that replace the path search by known algebraic paths and further reduce the computational complexity of the algorithm. Based on the symmetric subarrays, Zhi et al. [8] applied the generalized ESPRIT algorithm to search the azimuths of each near-field source. Liu Liang et al. [9] used the idea of rank reduction to construct a manifold matrix that was only related to the steering angle parameters and was used to search the near-field source azimuths. References [10][11][12] reconstructed a matrix related only to the azimuth parameters of the sources, and used the MUSIC algorithm to search for the azimuths. After obtaining the source azimuths estimation, the literatures [8][9][10][11][12] searched the range parameter in corresponding directions, which transformed the two-dimensional search into multiple one-dimensional searches, this type of method had a large array aperture loss. Some studies [13,14] estimated the source azimuths by constructing a higher order statistics (HOS) matrix to reduce the array aperture loss. There are also some other near-field parameter estimation methods, such as, the maximum likelihood estimator proposed in [15] and the weighted linear prediction method presented in [16].
With the successful application of sparse signal reconstruction in far-field DOA estimation [17] and the excellent performance of this kind of algorithm in anti-noise ability and snapshot number, more and more scholars also carry out the research on sparse reconstruction based near-field localization methods. Wang et al. [18] proposed a hybrid source localization algorithm based on the sparse representation of cumulants. Tian et al. [19] proposed a joint algorithm of MUSIC and sparse representation for localization of mixed sources with better estimation performance. In [20,21], the covariance matrix of the received signal of the virtual far-field array was constructed and the two-dimensional near-field parameter estimation problem was converted into the reconstruction problem of two one-dimensional sparse signals. Hu et al. [22,23] achieved a sparse estimation of DOA and range by sparse representation of anti-diagonal elements of received signal covariance matrix, which is similar to the method of [20], but with lower computational complexity. In [24,25], by using the sensor-angle distribution to characterize the sensor-dependent phase progression as a function of the source range and its direction, the sensor-dependent spatial frequency signature was estimated by sparse reconstruction techniques, and the results were then mapped back to source range and DOA estimation for the near-filed source localization.
The background noise of the above near-field localization algorithms assumed that it followed Gaussian distribution or at least had a finite second-order statistics, but in the natural environment or many engineering applications, the noise often exhibits non-Gaussian impulsive nature, and often has large amplitudes in a short time [26]. Mathematically, α-stable distribution model can be used to describe this type of noise [27]. Studies have shown that α-stable distribution process does not have the statistics of order greater than the characteristic exponent α(0 < α ≤ 2), so the algorithms based on the second-order statistics (SOS) and HOS will degrade or even fail in α-stable distributed noise environment. To this end, a number of articles had been studied, one way is to replace the SOS with fractional lower order statistics (FLOS), such as the fractional lower order correlation (FLOC) and the phased fractional lower order correlation (PFLOC) [28,29]. However, FLOS needs to know the characteristic exponent of the α-stable distribution in advance, which is difficult to obtain in practice. Another way is to replace the SOS with the robust second-order statistics. For example, the secondorder correlation was replaced by a robust correlation, nonlinear transform correlation (NTC) [30,31], which did not need to know the characteristic exponent of the α-stable distribution in advance. To solve the localization problem of near-field sources in the impulsive noise, Wang et al. [32] constructed matrixes by using FLOS and estimated position parameter of near-field sources by rooting method. By combing the concept of PFLOC and the GESPRIT method in [8], Qiu et al. [33] proposed a new search-free method for near-field source localization under impulsive noise, referred to as the search-free PFLOC-based GESPRIT. We studied the near-field localization problem based on sparse reconstruction of the constructed far-field approximate FLOC vector by exploring the symmetrical characteristic of the array under impulse noise environment for the first time in [34]. This paper is on this basis to do further research on this issue. In this paper, to avoid estimating the characteristic exponent α of the impulsive noise before using the FLOC statistics, we define a robust correlation vector, NTC vector, which is just like the FLOC vector that is in relation with the DOA only, then the bearings can be estimated based on the sparse reconstruction of the FLOC vector or NTC vector. With the estimated bearings, the range can consequently be obtained by the sparse reconstruction of the output of a virtual array. Lots of simulation experiments are made in this paper to demonstrate that the proposed algorithms have the merits of good noise suppression ability.

α-Stable Distribution
α-stable distribution is the only kind of distribution satisfying the generalized central limit theorem. Compared with the Gaussian distribution, the α-stable distribution has a thicker statistical tail, so it has significant impulse characteristics in time domain. There is no closed expression for the probability density function of the α-stable distribution, but the characteristic function can be used to facilitate the representation of the α-stable distribution It can be seen that the characteristic function of the α-stable distribution is determined by four parameters: α(0 < α ≤ 2) is the characteristic exponent, which describes the impulsive degree of the distribution. The smaller the α is, the thicker the corresponding distribution tail is, the more significant impulsive the signal is; β(−1 < β < +1) is the symmetry parameterand that β = 0 corresponds to the symmetric distribution, abbreviated as symmetry α-stable (SαS) distribution; γ(γ ≥ 0) is the dispersion which is similar to the variance of Gaussian distribution; a(−∞ < a < +∞) is location parameter, for the SαS distribution it represents the median or the mean. When α = 2 and β = 0, the characteristic function is the same as that of Gaussian distribution, that is, the Gaussian distribution is a special case of the α-stable distribution. A very important characteristic of the α-stable distribution is that it does not have a finite two order statistics and higher order statistics.

Signal Model
Consider the case of K independent narrowband sources s 1 (t), s 2 (t), · · · , s K (t) be in the near-field of a symmetric uniform linear array (ULA) with N = 2M + 1 isotropic sensors as illustrated in Figure 1. With the array center being the phase reference point, the signal received by the mth sensor at time t can be expressed as where n m (t) is the additive noise, and τ mk is phase shift related with the kth source signal's propagation time delay between phase reference point and sensor m. By Fresnel approximation, it can be given by w k and φ k with the following form Herein, λ is the wavelength, d is the interspacing,(θ k , r k ) represent the DOA and range parameters of the kth source.
The signal model in Equation (4) can be concisely expressed as where The purpose is to jointly estimate the DOAs θ 1 , θ 2 , · · · , θ K and range parameters r 1 , r 2 , · · · , r K for multiple near-field sources from the received array data X = [X(1), X(2), · · · , X(L)] where L is the number of snapshots. Assume the source signals are statistically mutually independent with zero mean and the noises n m (t) are complex isotropic Gaussian distributed random processes and are independent of the source signals, the SOCSR algorithm proposed in [23] localized the near-field sources by solving two spare reconstruction problems of the second order correlation vector of the array received signal. When the noises are complex isotropic SαS distributed random processes, the performance of the SOCSR algorithm will degrade since SαS distribution does not have finite SOS. This fact will be verified by the simulations in Section 5. In order to inhibit the influence of the α-stable impulsive noise, this paper proposes two new near-field sources localization algorithms by solving the spare reconstruction problem of the FLOC and NTC vector of the array received signal which is referenced as fraction lower order correlation-based sparse reconstruction method (FLOCSR) algorithm and nonlinear transform correlation-based sparse reconstruction method(NTCSR) algorithm.

FLOC Matrix of the Array Received Signal
The spatial fraction lower order cross correlation between the mth and nth sensor can be defined as [28,29] c FLOC (m, where In order to ensure c FLOC (m, n) is a finite value, the value of p needs to be less than the characteristic exponent α of the impulsive noise which is difficult to estimate in some practical applications.
The FLOC matrix of the array received signal can be formulated as where Λ = diag[Λ 11 , Λ 22 , · · · Λ KK ] and I is the identity matrix.

NTC Matrix of the Array Received Signal
In order to avoid estimating the characteristic exponent α of the impulsive noise in practical applications, a robust correlation, the nonlinear transform correlation (NTC), between the mth and nth sensor is defined as follows [30,31] where δ ≥ 1 is called scale factor. It has been proved that c NTC (m, n) is bounded and then the NTC matrix of the array received signal can be estimated.

Proposed Two Step Estimation Method
The FLOC c FLOC (m, n) and NTC c NTC (m, n) can be unified written as the robust correlation c x (m, n) and the corresponding matrix form can be written as C. When m = −n the robust correlation c x (m, n) is independent of the parameter φ k . This means that by exploiting the robust correlation between the symmetric sensors, we can transform the original two-dimensional (DOA and range) estimation problem into a one-dimensional (DOA) estimation problem. Stacking Equation (9) or Equation (11) for the symmetric sensors, we can build a virtual far-field model where −1), . . . , c x (M, −M)] T ∈ C 2M×1 and Λ w = [Λ 11 , Λ 22 , · · · , Λ KK ] T is the received signal vector and source signal vector of the virtual far-field array, the manifold matrix of the virtual far-field array can be expressed as A w (θ) = [a w (θ 1 ), a w (θ 2 ), · · · , a w (θ K )] ∈ C 2M×K with the virtual array steering vector a w (θ k ) = e −j2Mw k , · · · , e −j2w k , e j2w k · · · , e j2Mw k T ∈ C 2M×1 .

Step-1: DOA Estimation
The virtual far-field array received signal vector c x can be sparsely represented in a redundant basis. Define a setθ = θ 1 ,θ 2 , · · · ,θ N θ which denotes potential DOAs of interest sources and assume that the true DOAs are exactly on this set. The number of the potential DOAs N θ should be much greater than K which is the number of sources. Define the overcomplete basis A w (θ) = a w (θ 1 ), a w (θ 2 ), · · · , a w (θ N θ ) and the potential source signal vector v = [v 1 , v 2 , · · · , v N θ ] T . As a result c x can be rewritten as the following form It can be seen that the elements of vector v have K nonzeros, that is v k = Λ ii ifθ k = θ i , i = 1, · · · , K. Hence the DOA estimation problem can be reduced to finding the nonzero elements of the vector v.
Since v is sparse, so it can be estimated by solving the following sparse reconstruction problem where ε 1 is a parameter which means how much of the error we wish to allow and plays an important role in the algorithm performance.

Step-2: Range Estimation
Given a DOA θ k (k = 1, · · · , K) which is estimated in Step1, the robust correlation matrix of the near-filed received signal can be written as where the manifold matrix A θ k , r = a θ k , r 1 , a θ k , r 2 , · · · , a θ k , r K . Applying the vectorization operator on Equation (15), we have B θ k (r) = a * θ k , r 1 ⊗ a θ k , r 1 , · · · , a * θ k , r k ⊗ a θ k , r k ∈ C N 2 ×K where ⊗ denotes Kronecker product. It is interesting to see that y θ k in Equation (16) can also be regarded as output of the virtual far-filed array where B θ k (r), Λ w and ξvect(I) are the virtual manifold matrix, equivalent source signal vector, and equivalent noise vector, respectively. Notice that the vector ξvect(I) has only N nonzero elements, then these elements of y θ k corresponding to these positions of nonzero elements in ξvect(I) can be removed and the rest N(N − 1) entries of y θ k corresponding to these positions of zeros elements in ξvect(I) can be preserved. Then, the virtual far-field array output vector processed by this operation can be formulated as where B θ k (r) ∈ C N(N−1)×K is the new manifold matrix obtained by removing the rows of matrix B θ k (r) which are corresponding to these positions of nonzero elements in ξvect(I). This elimination operation can further reduce the effect of the impulsive noise. Using a similar approach as in Step-1, the virtual received signal vector y θ k can be sparsely represented as the following form y where B θ k (r) ∈ C N(N−1)×N r is the overcomplete basis on a setr = [r 1 ,r 2 , · · · ,r N r ]. N r is the number of the potential sources on the direction of θ k and should be much greater than the number of real sources N r k on the direction of θ k , p = [p 1 , p 2 , · · · , p N r ] T is the potential source signal vector that have N r k nonzeros, that is p j = Λ k i k i ifr j = r k i i = 1, · · · , N r k . Hence the range parameter estimation problem can be resolved by finding the nonzero elements of vector p, which can be estimated by solving the sparse reconstruction problem given by where ε 2 is also a parameter which means how much of the error we wish to allow and plays an important role in the algorithm performance.

Simulation Results
In order to verify that our proposed FLOCSR and NTCSR algorithms have a more noise-rejection ability than the SOCSR algorithm in the α-stable impulsive noise environment, we conduct a series of numerical experiments under a variety of simulation conditions. In order to solve the convex optimization problem given in Equations (14) and (20), the software package CVX [35] is used. Considering an N = 15 element ULA, the separation distance between the elements is λ/4. Uniform sampled in the angular space [−90 o , 90 o ] at 1 o interval and near-field range scope [0.1λ, 10λ] at 0.1λ interval, that is N θ = 181 and N r = 100.
As the characteristic of the α-stabledistribution makes the use of the standard SNR meaningless, a new SNR measure, generalized signal-to-noise ratio (GSNR) is defined as [27] GSNR = 10 log 10 σ 2 where σ s is the variance of the signal. Two performance criteria are used to assess the performance of the algorithms. The first one is the probability of success. A successful simulation is defined if angle difference between the estimated DOA and the real DOA is less than 3 o and distance difference between the estimated range and the real range is less than 0.3λ for all incident sources. The probability of success is defined as the ratio of the number of successful simulations to the total number of Monte Carlo simulations. Another criterion that used to assess the performance of the algorithms is the average root mean square error (RMSE) defined as where x l is the real value of the DOA or the range parameter, x l (i) is the ith estimation value of x l and N ok is the number of successful simulations.

Simulation 1
Two independent Gaussian distribution sources with equal power in the near-field region at locations (θ 1 , r 1 ) = (20 o , 1.5λ) and (θ 2 , r 2 ) = (45 o , 3.6λ) are considered. Suppose the noise is modeled to be SαS distributed with α = 1.5. The generalized SNR is set to be GSNR = 5 dB and the number of snapshots is L = 500. Two-hundred Monte-Carlo simulations are performed individually for each method. Figures 2-4 give the estimated locations of the SOCSR, FLOCSR, and NTCSR algorithm. It can be seen that the simulation results of the SOCSR algorithm are scattered around the real location, whereas the simulation results of the FLOCSR and NTCSR algorithms are closely gathered near the real location, especially the NTCSR algorithm.

Simulation 2
In this simulation the locations of two independent Gaussian distribution sources with equal power are set as (θ 1 , r 1 ) = (20 o , 0.2λ) and (θ 2 , r 2 ) = (45 o , 1.4λ). The characteristic exponent of the α-stable noise is fixed at α = 1.5 and thenumber of snapshots is L = 1000. Figure 5 shows the performance of the three algorithms for various GSNRs ranging from 2 dB to 20 dB. We see that probability of success of all algorithms improve with the increase of GSNR, and the proposed FLOCSR and NTCSR algorithms outperform the SOCSR algorithm. For example, when GSNR = 6 dB, the probability of success of SOCSR algorithm is only 58% whereas the probability of success of FLOCSR and NTCSR method are all above 90%. The RMSE of DOA and range of all algorithms decrease with the increase of GSNR but the proposed FLOCSR and NTCSR algorithms have a less value than the SOCSR algorithm at the same GSNR. That is to say, the proposed FLOCSR and NTCSR methods have a better estimation accuracy and precision than SOCSR algorithm. From the simulation results, we can also see that the performance of the NTCSR algorithm is better than that of the FLOCSR algorithm, indicating that the nonlinear transform correlation has a better ability to suppress impulse noise compared with the fraction lower order correlation.   Figure 6 plots the performance of the three algorithms varying with different values of the characteristic exponent of the α-stable impulsive noise. The simulation environment is as same as Simulation 2 except that the GSNR is kept at 10 dB. As shown in Figure 6, our proposed FLOCSR and NTCSR algorithms demonstrate their performance enhancement over SOCSR algorithm in the sense of the probability of success and RMSE under the highly impulsive noise environment. In particular, the performance of NTCSR algorithm which has a good ability to suppress impulse noise is more outstanding when in the strong impulse noise environment with characteristic exponent α = 1.1, the probability of success is almost 1 and the lowest RMSE.

Simulation 4
In this experiment, the simulation environment is accordance with Simulation 2 except for GSNR = 6 dB and the number of snapshots is varied from 50 to 1600. Figure 7 shows the relationship between the performance of the three algorithms and the number of snapshots. We can observe that as the number of snapshots increase all algorithms exhibit a decrease in RMSE and an increase in the probability of success. However when the number of snapshots is more than 400, the increased performance caused by the snapshots is not obvious. Nevertheless, our proposed FLOCSR and NTCSR algorithms produce lower RMSE and higher probability of success compared to the SOCSR algorithm when the same number of snapshots is used.

Simulation 5
The ranges of two incident sources are fixed at r 1 = 0.2λ, r 2 = 1.4λ, Figure 8 shows the capability of angular separation of three algorithms with the DOA of the first source is fixed at θ 1 = 20 • and the DOA of the second source is varied from 24 • to 40 • with an interval of 2 • under the simulation environment of GSNR = 6 dB, α = 1.5, and L = 1000. It is generally considered that the greater the angle separation, the smaller the influencesbetween the two incident sources, the better the estimation performance. The simulation results verify this point of view, that the greater the angle separation, the higher the probability of success and the lower of the RMSE of the angle estimation. Since the range parameter is obtained on the basis of the estimated DOA, thence, the performance of the range estimation has little to do with the angle separation. Under the same angle separation condition, since the SOCSR algorithm is not resilient to impulsive noise, it realizes the localization at a lower probability of success and a higher RMSE relative to the FLOCSR and NTCSR algorithm which can effectively suppress the impulse noise.

Simulation 6
To test the capability of range separation of the algorithms, we fix the directions of two incident sources at θ 1 = 20 • , θ 2 = 40 • and the range of the first source at r 1 = 0.2λ. The range of the second source is varied from r 2 = 0.2λ to r 2 = 2λ with an interval of 0.2λ. The plots in Figure 9 are obtained under GSNR = 6 dB, α = 1.5 and L = 1000. According to the principle of the algorithm, DOA estimation is independent of range, accordingly, the range separation changes will not affect the estimation results of DOA. Therefore, as shown in Figure 9, the probability of success and the RMSE of DOA hardly vary with range separation. However, the RMSE of range improves with the increase of range separation. In any case, our proposed two algorithms have better probability of success and lower RMSE than the SOCSR algorithm.

Simulation 7
In the NTCSR algorithm, the scale factor δ determines the degree of nonlinear transformation of the array received signals, in other words, it determines the degree of the impulse noise suppression. Figure 10 shows the performance of NTCSR algorithm in α = 1.2 and 1.5 two different impulse noise conditions for various scale factors ranging from 1 to 10. The number of snapshots is L = 1000 and GSNR = 10 dB. Although the localization of the near-field signal sources can be relatively successfully implemented in both impulse noise environments, different scale factors result in different estimation accuracies. From Figure 10 we can observe that δ ∈ [3 − 5] would be the optimal domain for NTCSR algorithm to achieve its best performance in the RMSE of DOA and range.

Conclusions
In this paper, the localization problem of near-field sources under impulsive noise environments is studied. Based on the symmetrical characteristic of the array, we construct the robust far-field approximate correlation vector in relation with the DOA only, which allows for bearing estimation based on the sparse reconstruction of the robust correlation vector. With the estimated bearing, the range can consequently be obtained by the sparse reconstruction of the output of a virtual array. Simulation results indicate the superiority of the presented two algorithms in the probability of success and RMSE under a variety of impulsive noise environments.