Structure-Aware Bayesian Compressive Sensing for Near-Field Source Localization Based on Sensor-Angle Distributions

A novel technique for localization of narrowband near-field sources is presented. The technique utilizes the sensor-angle distribution (SAD) that treats the source range and direction-of-arrival (DOA) information as sensor-dependent phase progression. The SAD draws parallel to quadratic time-frequency distributions and, as such, is able to reveal the changes in the spatial frequency over sensor positions. For a moderate source range, the SAD signature is of a polynomial shape, thus simplifying the parameter estimation. Both uniform and sparse linear arrays are considered in this work. To exploit the sparsity and continuity of the SAD signature in the joint space and spatial frequency domain, a modified Bayesian compressive sensing algorithm is exploited to estimate the SAD signature. In this method, a spike-and-slab prior is used to statistically encourage sparsity of the SAD across each segmented SAD region, and a patterned prior is imposed to enforce the continuous structure of the SAD.The results are then mapped back to source range and DOA estimation for source localization. The effectiveness of the proposed technique is verified using simulation results with uniform and sparse linear arrays where the array sensors are located on a grid but with consecutive and missing positions.


Introduction
Localization of emitters is an important topic with many applications in radar, sonar, and wireless communications.The majority of research work in array signal processing has been focused on far-field signals; that is, the source signals are assumed to be located far from the receiver.Accordingly, the propagation waves assume a planar wavefront, making source localization solely based on the direction-of-arrival (DOA) of the sources.However, this characterization becomes invalid when sources are located in the near-field region.In this case, each source is characterized by both DOA and range that vary with sensor positions.To solve near-field source localization problems, subspace-based techniques, for example, MUSIC, were extended to a two-dimensional field [1].However, this entails a high computational load required for both DOA and range searching.
To decouple these two variables of the DOA and range for simplified computations, a more accurate quadratic approximation of the wavefront, in lieu of the linear approximation, can be made with respect to the sensor positions [2].Consequently, the phase delay is no longer linear with the sensor positions but is decided by the DOA and range of the sources jointly.A number of methods have been proposed to jointly estimate the DOA and range information for nearfield source localization.For instance, algorithms based on high-order statistics were proposed in [2][3][4].The work in [2] exploits the rotational invariance in certain cumulantdomain signal subspace for range and bearing estimation.The technique in [3] applies a symmetric uniform linear array (ULA) and uses eigenvalues together with the corresponding eigenvectors of two properly designed matrices to jointly estimate signal parameters.As such, it avoids any spectral peak searching.Two orthogonally polarized sensors were utilized in [4] to improve performance.These methods, however, do not meet the desired performance and require certain array properties, for example, a ULA with an interelement spacing less than half a wavelength.In practical applications, 2 International Journal of Antennas and Propagation such approaches are very costly due to receiver hardware and computational complexity.
Sparse arrays [5][6][7] use fewer elements to achieve the same array aperture of fully populated arrays.They have similar mainbeam properties and, therefore, provide similar performance in terms of angular accuracy, resolution, and detection of targets close to interference directions, all being achieved with reduced size, weight, power consumption, and cost.An algorithm was proposed in [8] to cope with the near-field source localization problem using sparse arrays.However, it relies on the assumption that the sparse array configuration must be symmetric.In addition, similar to [2][3][4], the second-order Taylor approximation of the range between the source signals and the sensors is assumed to decouple the two variables.As a result, the performance is compromised when the sources are placed close to the array.
In this paper, we propose an alternative approach using sparse arrays by mapping the source range and DOA information into sensor-dependent phase progression.The sensor-angle distribution (SAD) is used to represent the phase progression in the joint space and spatial frequency domain with a signature that resembles those of instantaneously narrowband frequency-modulated (FM) signals in the joint time-frequency domain [9,10].In this respect, the signal processing algorithms introduced for quadratic timefrequency analysis can be readily applied for the processing of the SAD for near-field source localization.It is important to note that, while the SAD was originally developed in the form of spatial Wigner distribution [9,10], which is effective only for second-order (linear FM) spatial frequency signatures, we consider it in the general form within Cohen's class with a proper time-frequency kernel for improved sensor-angle relationship [11].The proposed technique permits the sensordependent phase progression as a high-order polynomial phase signal (PPS) to account for close sources, whereas the coefficients are estimated from the SAD signature with simple least squares fitting, which does not require parameter searching.
Another important contribution of this paper is the consideration of the source localization problem in a sparse array configuration where the array sensors are located in a uniform linear grid but with missing positions.As such, conventional time-frequency analysis tools become infeasible due to the undesirable artifacts caused from missing sensor entries [12].Instead, we use the emerging compressive sensing (CS) [13] techniques to reconstruct the SAD in the joint space and spatial frequency domain.The problem is considered in the Bayesian compressive sensing (BCS) framework that generally achieves desired performance and enables flexible utilization of various signature structures through the application of proper priors [14][15][16][17].In particular, we use the structure-aware BCS, which was developed for effective reconstruction of sparse time-frequency distributions of FM signals, to exploit the sparsity and continuity of the spatial frequency signatures in the SAD domain.In doing so, the isolated or sporadic entries in the reconstructed SAD caused by the presence of missing sensors and measurement noise can be significantly suppressed, leading to more reliable SAD reconstruction and spatial frequency estimation.The obtained results are then mapped back to the estimation of source range and DOA information for effective source localization.Compared to the approach proposed in [8] using a sparse symmetric array, the proposed technique offers higher flexibility and robustness to array design and improves the source localization accuracy when the sources are close to the array and thus the phase exhibits a high-order polynomial phase relationship with the sensor positions.
The rest of the paper is organized as follows.Section 2 introduces the signal model of near-field localization based on a third-order approximation.The SAD and its sparse reconstruction are described in Section 3, and the structureaware Bayesian compressive sensing technique applied for sparse SAD reconstruction and spatial frequency estimation is described in Section 4. Simulation results are provided in Section 5 to numerically evaluate and compare the performance.Such results reaffirm and demonstrate the effectiveness of the proposed technique.Section 6 concludes this paper.
Notations.We use lower-case (upper-case) bold characters to denote vectors (matrices).In particular, I  denotes the × identity matrix.(⋅) * denotes complex conjugate, whereas (⋅)  and (⋅)  , respectively, denote the transpose and conjugate transpose of a matrix or vector.diag(x) represents a diagonal matrix that uses the elements of x as its diagonal elements.‖ ⋅ ‖ 2  2 implies the Euclidean ( 2 ) norm of a vector.  (⋅) expresses the probability density function (pdf) and ( | −) denotes the conditional pdf of random variable , given other parameters.Beta(, ) computes the beta function for corresponding elements  and .CN( | , ) describes that random variable  follows a complex Gaussian distribution with mean  and variance .Furthermore, () is the Dirac delta function of .

Signal Model
For the simplicity of presentation, the signal model is presented for the single-source case, but the proposed technique can be straightforwardly extended to the multiple source problems and simulation results are provided for two source scenarios.
Consider a near-field narrowband source impinging on a uniform or sparse linear array with  = 2 + 1 elements, as depicted in Figure 1, where  ∈ Z + .Note that we assume  to be an odd integer without loss of generality, but the problem can be similarly formulated for an even value of .Examples of even values of  are provided in Section 5. Let the index of the reference element be 0, and the sensors are located at   , where   ∈ Z for  = −, . . ., , and  is the unit interelement spacing.To avoid spatial ambiguity,  is usually taken as half-wavelength, denoted as  = /2.The source is located with a range  and azimuth angle  with respect to a reference sensor.Note that the array may be asymmetric; that is, it is not necessary to assume  − = −  , for  = −, . . ., .In a near-field environment, the range between the source and the sensors varies for each of the  sensors.It is straightforward to show that the range between the source signal and sensor  is given by

Sensors
As such, the near-field steering vector is given by For source signal (), the data received at the sensor array is expressed as where n() is the additive Gaussian noise vector and  is the number of observed time samples.
To understand the significance of the steering vector, consider the special case of fully populated array that is uniformly spaced with interelement spacing ; that is,   = .Then, (1) can be expressed as Assume that  is reasonably larger than .Then, using the Taylor series expansion, we can obtain the third-order approximation, which is given by It is indicated that the near-field steering vector given in ( 2) is a uniformly sampled FM signal with respect to .
In this respect, the signal processing techniques devised for nonstationary signal analysis and quadratic time-frequency signal representations can be readily applied for near-field source localization in the joint space and spatial frequency domain.The corresponding spatial frequency is obtained from the phase difference with respect to , given as where It is clear that the polynomial coefficients  0 ,  1 , and  2 are related to the range and the DOA.As such, these parameters can be estimated once the sensor-dependent spatial frequency is obtained.Note that, for a large range ( ≫ 1), () behaves as a chirp when the second-order term in ( 6) is ignored.Furthermore, when the source is located in the far field, as  approaches infinity, () approaches a constant, implying a single spatial frequency associated with a sensor-independent DOA.According to (7), from the estimated values of  0 ,  1 , and  2 , the unknown range and DOA can be, respectively, estimated as Specifically, when  = /2, the above expressions become

Sensor-Angle Distribution and Sparse Reconstruction
Several methods are available for the estimation of linear FM and PPS parameters (see, e.g., [18][19][20]).Such methods, however, are effective only when the samples are uniformly applied (i.e., the array is uniform linear in the underlying problem).When the array is sparse, as considered in this paper for more effective array aperture utilization, these methods become difficult to apply.For illustrative purposes, we consider a sparse linear array structure, consisting of a set of uniform linear subarrays with different coprime interelement spacings, as an example.The coprime array structure was proposed in [21] for systematical sparse array configuration, which generally achieves ( 2 /2) coarray sensors with  physical sensors.Therefore, the aperture and International Journal of Antennas and Propagation the number of degrees-of-freedom (DOFs) can be improved using correlation-aware techniques.Coprime arrays have found broad applications in DOA estimation and interference suppression [22][23][24][25].
For the spatial frequency estimation, we exploit the SAD to describe the relationship between spatial frequency and the sensor positions.The SAD draws parallel to the timefrequency distribution.In essence, the SAD uses sensor positions in lieu of time samples whereas the spatial frequency replaces temporal frequency.With this analogy, SAD becomes a tool to reveal the local behavior of the spatial frequency along the space axis with respect to the array sensor positions.For sparse arrays, the SAD suffers from missing artifacts due to missing samples.We use time-frequency kernel and sparse reconstruction to achieve satisfactory SAD reconstruction and spatial frequency estimation.In the following, we describe the SAD definition in Section 3.1 and discuss the effect of array sparsity in Section 3.2.The sparse reconstruction of SAD from a sparse array is addressed in Section 3.3, and the structure-aware Bayesian compressive sensing that implements the sparse reconstruction is detailed in Section 4.

Sensor-Angle Distribution.
The SAD corresponding to the th sensor position and spatial frequency  at time instant  is expressed as where (; ) denotes the th element of vector x() and (V, ) is a kernel function that characterizes the distribution and is a function of sensor position and sensor lag [26].The above equation defines Cohen's class of power distribution over the joint variables  and .In this respect, all the standard kernel designs applied in the time-frequency analysis literature can be employed with the SAD.The distribution in (10) can be averaged over the  samples to reduce the effect of observation noise, yielding Stacking (, ) with respect to  and  forms an SAD matrix D.

Data-Dependent Time-Frequency Kernel.
Benefiting from the analogy of SAD to time-frequency distributions, we understand that the sparse placement of array sensors is deemed to produce artifacts that clutter the SAD, making the spatial frequency estimation difficult.Analysis of the effects of missing samples on time-frequency distribution is given in [12].To mitigate effects of the induced artifacts and improve time-frequency signature estimation, the first step is to apply a time-frequency kernel [12,17,27].Time-frequency kernel allows mitigation of cross-terms between different signal components as well as the effect of artifacts due to missing samples (dual of missing sensors in the underlying problem).Time-frequency kernels are generally represented as a multiplying weight (mask) matrix in the ambiguity domain.The ambiguity function corresponding to the underlying SAD is defined with respect to sensor position lag  and spatial frequency shift  as The averaged AF is then obtained by averaging over the  time samples as The time-frequency kernels are based on the fact that majority of autoterm distributions are located near the origin in the ambiguity domain, whereas the cross-terms and artifacts are distant from the time lag and frequencyshift (sensor lag and spatial frequency shift in the underlying problem) axes.This property has motivated researchers to propose educed-interference distribution kernels with lowpass filter characteristics to suppress cross-terms and preserve autoterms.Such kernels can be signal-independent or signal-dependent.In this paper, we use the adaptive optimal kernel (AOK) and modify it to specifically deal with missing positions.The AOK is formulated as an optimization problem in [28]: where C denotes a matrix with all (, ) entries.Note that the ambiguity function (, ) and the kernel (, ) are defined in the polar coordinates with the radius  and angle .The objective of the AOK is to maximize the energy in the ambiguity domain in order to preserve the autoterms.At the same time, the Gaussian constraint in (, ) leads to cross-term reduction.The choice of  influences the tradeoff between cross-term suppression and autoterm preservation.

Sparse Reconstruction of SAD from Sparse Array.
As addressed above, a key step to achieve improved SAD reconstruction and spatial frequency estimation in the presence of missing sensors is through sparse reconstruction.In extending the techniques developed for sparse reconstruction of time-frequency distribution in [12,17,27] to the underlying sparse SAD reconstruction, we define the following discrete sensor autocorrelation function (SAF) for time instant : The averaged SAF can be similarly defined using all available samples to reduce the effect of additive noise.The SAF matrix A is formed that contains (, ) for all the entries of  and .As such, at each sensor position, the columns of the averaged SAD matrix D and that of the averaged SAF matrix A are related by a onedimensional (1D) Fourier transform along the lag dimension.
Denote a [] as a column vector that contains all the elements in the th column of A and d [] as a vector corresponding to the th column of D. Then these two vectors are related by a [] = Fd [] ,  = 1, . . ., , where F is the inverse Fourier transform matrix.Equation ( 17) describes a linear model relating the sparse vector d []  to the compressed observation vector a [] using the inverse Fourier dictionary.As such, d [] can be obtained by solving a sparse optimization problem, which is commonly used in compressive sensing techniques.The positions of nonzero entries in d [] indicate the estimations of spatial frequency.Note that F can be designed as a wide matrix so as to achieve a higher spatial frequency resolution in the estimated d [] .As demonstrated for time-frequency analysis, the 1D relationship between the instantaneous autocorrelation function (IAF, the dual of SAF) and the time-frequency distribution (the dual of SAD) is much more effective than the conventional two-dimensional (2D) approach between the ambiguity function and the time-frequency distribution [29,30] and allows local time-frequency reconstruction for each time instant (the dual of sensor positions).In particular, such local reconstruction capability is fundamental to enable the development of the structure-aware Bayesian compressive sensing technique to be described in Section 4.

Structure-Aware Bayesian
Compressive Sensing

Continuous Structure Modeling with Spike-and-Slab Prior.
A number of approaches, such as the OMP and LASSO [31,32], can be effectively applied for the reconstruction of sparse SAD signatures.To exploit the continuous structure of the SAD signature, however, the BCS-based approaches are considered most effective due to their flexibility of employing proper priors for this purpose.We use the structure-aware BCS technique [17] to improve estimation performance.
To encourage sparsity of the SAD signatures, we place a spike-and-slab prior to each element of d []  [33]; that is, for its th element, where  []   is the prior probability of a nonzero element and  []   is the precision (reciprocal of variance) of the complex Gaussian distribution, for the th sensor position.
To facilitate an analytical inference, which involves the delta function, the product of following two latent variables is introduced to follow the pdf in (18): where is a binary variable following the Bernoulli distribution [17].
In addition, a Gaussian prior is placed on the additive noise as  []   ∼ CN( []   | 0,  −1 ), where the precision  is assumed to be independent of  and .To acquire a trackable posterior distribution, we place Gamma priors, which are conjugate to the Gaussian distribution, on both  []   and .To encourage continuity in the joint space and spatial frequency domain, we utilize data in the neighborhood sensor positions between  − 1 and  + 1, when the SAD signature in sensor position  is estimated.As shown in Figure 2, we categorize the relationship into three different patterns, termed as Pattern 1 ("strong rejection"), Pattern 2 ("weak rejection"), and Pattern 3 ("strong acceptance"), respectively.Define z1 , . . ., z6 as illustrated in Figure 2. Based on the idea of continuity, we define two variables,  []    = ∑ 6 =1 z and  []   = ∑ 3 =1 z z7− , which, respectively, represent the number of nonzero neighboring entries and the number of nonzero diagonal pairs.It should be noted that the vertical pixel pairs are excluded for being impractical in the underlying problem.
Pattern 1 includes two cases of  []    = 0 or  []    = 3.In the former case, all neighbors are zero values, and thus the th entry would be zero valued with a high probability.In the latter cases, all the neighbors are nonzero values.In both cases,  0 <  0 is assumed in prior Beta( 0 , 0 ) to encourage a small value of  []   so as to reject this entry.Pattern 2 illustrates those cases of  []    > 0 or  []    ̸ = 1. Figure 2(b) shows one of the examples with a single nonzero neighboring entry.In this case, the probability that the th entry takes zero values is fair, and  1 =  1 is used in the prior Beta( 1 , 1 ) to exert noninformative prior on  []   .
International Journal of Antennas and Propagation shows two examples of such cases.We confidently believe that the th entry would take a nonzero value with a high probability when nonzero neighboring entries have such continuous and symmetric structure in the 2D space and spatial frequency domain.In this case,  2 >  2 in the prior Beta( 2 , 2 ) is assumed to encourage a large value of  []   to accept this entry.
Note that the above patterns encourage instantaneously narrowband spatial frequency signatures and closely follow those developed for time-frequency signature reconstruction of FM signals [17].These patterns are, however, different from those used to encourage merely spatial connectivity as described in [34,35].

Bayesian Inference.
The Bayesian inference of the proposed algorithm is carried out by the Gibbs sampler [33].The paired Gibbs sampler iteratively samples the following conditional pdf:  ( []   ,  []   |  []  \ , z []  \ , a [] ) =  ( []   |  []   ,  []  \ , z []  \ , a [] ) ×  ( []   |  []  \ , z []  \ , a [] ) , where  []  \ and z []  \ , respectively, denote  []   except the variable  []   and z [] except the variable  []   .The probability ( []    = 1 |  []  \ , z []  \ , a [] ) is acquired analytically as where  []   is derived as with , and F  represents the th column of the inverse Fourier transform matrix F. The conditional distribution of ( []    |  []    = 1,  []  \ , z []  \ , a [] ) is expressed as For  []    = 0, the value of  []   does not affect the result of  []   .Thus, we can draw the value of variable  []   from its prior.We adopt the procedure in [17] to update the mixing weight  []   , signal precision  []   , and noise precision , in order to obtain  []   and  []   .Then, the sparse vector d [] , whose positions of nonzero entries indicate the estimated nonzero spatial frequencies () in f [] , can be obtained using (19).
Based on (6), the following linear model can be obtained: where , . . ., ()  ]  , and As such, once the sparse spatial frequencies along the array axis are estimated, the polynomial phase coefficients in b can be obtained using least squares fitting as Then, the range and DOA can be estimated based on (8) or (9).

Simulation Results
Simulation examples are presented to demonstrate the effectiveness of the proposed technique.For each sensor, the maximum number of iterations in the Gibbs sampling is 200, and the sampler with the maximum marginal likelihood in the last 20 samples is chosen to estimate d [] for  = −, . . ., .
For all examples, the carrier frequency is considered to be 15 MHz, yielding a wavelength of  = 20 m.As such, the unit interelement spacing is set to  = /2 = 10 m.The source signals are assumed to be quadrature phase-shift keying (QPSK) modulated.In all simulations, we use a 32-element ULA as well as its subset described as a linear sparse array with missing sensors.For the latter, as shown in Figure 3,  It should be noted that this configuration is not symmetric and, as such, the approach proposed in [8] becomes inapplicable.We evaluate the source localization performance through Monte Carlo simulations.The average root mean square error (RMSE) of the estimated source location is used as the performance metrics, expressed as where   order to obtain a smooth figure, we oversample the spatial frequency on a multiple of 128. Figure 4 shows the SAD results for a ULA from both Wigner-Ville distribution (WVD) and sparse reconstruction, and the corresponding results for the sparse array are illustrated in Figure 5. Figure 4(a) shows the WVD of the SAD.Because there is no missing sensor, the signature can be estimated using a number of conventional methods as well as sparse reconstruction techniques, with respect to the SAF.The corresponding results obtained from the OMP and the structureaware BCS are depicted in Figure 4(b).It is apparent that the signature is well approximated as a chirp and there is little difference noticed between the quadratic and thirdorder approximations.As such, in this case, the use of  The corresponding 16-element sparse array scenario is presented in Figure 5, in which the conventional method developed in [8] is not effective.In addition, the SAD is contaminated with artifacts that may obscure the true feature estimation due to missing sensors, as shown in Figure 5 It is evident that a desired performance is achieved.Figure 6 compares the RMSE performance as a function of the input SNR obtained from 100 independent trials.It is evident that the source localization performance is improved as the input SNR increases.The WVD in the ULA case presents the best performance since the WVD is known to be optimal in the analysis of linear FM signals [36,37].However, the performance of WVD degrades when the spatial frequency is nonlinear to the sensors positions or when some sensors are missing.The sparse array with the proposed method yields a similar performance but suffers from a gap of performance degradation.Such performance degradation partly stems from the signal power loss with the fewer number of sensors.

Example II:
Single-Source with a Close Range.In this example, the source is located at a range of 350 m and DOA of 30 ∘ with respect to the reference sensor.The corresponding coordinate is [303.1,175] m.As a consequence, the signature behaves as a PPS with respect to the sensor positions.The WVD and the estimated results obtained from the ULA and sparse array are shown in Figures 7 and 8, respectively.As indicated in Figure 7, the estimated spatial frequency is approximated as a third-order PPS with respect to the sensor positions.In this case, the WVD exhibits some cross-terms.The third-order PPS model yields more accurate results as compared to the quadratic approximation in this case to reflect the closer source range to the array.Under the third-order PPS approximation, the coefficient b = [0.4277,−0.004, −0.0001]  is obtained, resulting in the source position [286.14, 173.16] m, a 17.06 m shift from the true source position.For the sparse array, the effect of artifacts is again confirmed in Figures 8(a) and 8(b).This effect is significantly suppressed after the AOK is applied, as shown in Figures 8(c) and 8(d).
In Figure 9, we compare the RMSE performance of the second example.It is shown that the AOK-based results achieve better performance because of the effective crossterm and artifact suppression.In addition, there appears a floor in the obtained RMSE performance due to the spatial frequency resolution and bias.For the sparse array, a similar performance is attained to that of the ULA when the input SNR is higher than 0 dB and outperforms the WVD-based ULA results.

Example III:
Two-Source Scenario.Finally, we consider a two-source case in Figure 10.The first source is located at  = 300 m and  = 30 ∘ with respect to the reference sensor, whereas the second source is located at  = 1000 m and  = 60 ∘ .For the 32-element ULA, as shown in Figure 10(a), the SAD in this case is contaminated by the cross-terms between the two sources, making the spatial frequency estimation difficult.Such cross-term effect can be substantially reduced by exploiting the AOK, as shown in Figure 10(b).The sparse reconstruction results of OMP obtained from the AOK are depicted in Figure 10(c), where the spatial frequencies are not correctly identified towards both edges.Improved results are obtained from the proposed technique exploiting the structure-aware BCS, as depicted in Figure 10(d), due to the enforced spatial frequency continuity in the sparse

Figure 1 :
Figure 1: The geometry for near-field source localization.

Figure 2 := 1 .
Figure 2: Three patterns for 2D image.White squares denote entries with zeros amplitudes, blue shaded squares denote the entry with nonzero amplitudes, and red shaded squares denote the entry under test.(a) Pattern 1: strong rejection; (b) Pattern 2: examples for weak rejection with  []  = 1 case; (c) Pattern 3: examples for strong acceptance with  []  = 1 case.