Target Localization Based on Distributed Array Networks With Magnitude-Only Measurements

The source localization problem based on distributed array networks is formulated into a group sparsity-based phase retrieval problem where only the magnitude of received signals is available. Under such a framework, a 2-D localization method is proposed, where unlike traditional methods, random phase errors at array sensors will not affect estimation results. In addition, to deal with the off-grid problem for sparsity-based solutions, a model for off-grid bias is proposed and an efficient two-step method is developed accordingly to solve the 2-D off-grid problem. Simulation results show that the proposed solutions can solve the problem effectively.

For AOA-based methods, a structure with multiple sensor arrays distributed in a 2-D space is employed, where synchronization among all distributed sensor arrays is not required.There are normally two steps: the first is applying existing direction of arrival (DOA) estimation methods, such as those proposed in [11], [12], and [13], to estimate AOAs at all distributed sensor arrays, whereas the second is to find intersections of those estimated AOAs in order to localize the sources [9], [10].Recently, in [14], with the distributed array network, information across all sensor arrays is jointly exploited and the source localization problem was reformulated into a group sparsity maximization problem, where the area of interest is divided into grids along the x-axis and y-axis; a common spatial sparsity support corresponding to all distributed sensor arrays is enforced, leading to a better estimation performance, which also avoids the possible pairing and ambiguity problems associated with a two-step AOA-based solution.
The abovementioned AOA-and sparsity-based methods have assumed that there are no phase errors in the array model.In real applications, however, the phase information may not be reliable due to various reasons, and in the extreme case, the phase information may be lost completely, which unavoidably leads to an inaccurate estimation result.On the other hand, inspired by the work in noncoherent (i.e., magnitude-only measurements) DOA estimation, where only magnitude information is captured at all sensors [15], [16], [17], [18], [19], [20], [21], [22], a noncoherent source localization problem has been proposed in our earlier conference publication [23].The ambiguities associated with the noncoherent measurements at each subarray are resolved by employing a uniform circular array (UCA) [24]; otherwise, some reference source signals may be needed to remove the resultant ambiguities [16], [21].
One problem associated with the sparsity-based estimation approach is that the source locations are implicitly assumed to fall on the discrete grid points.However, in practice, quite often the assumption would not be satisfied, which leads to an off-grid problem.A straightforward solution is to employ a sufficiently dense grid so that the assumption can be considered to be valid to a great degree.However, this will increase the complexity significantly.Another solution is grid refinement, where instead of creating a dense grid initially to alleviate the off-grid problem, a coarse grid is first made, and then a denser grid is built around the estimated locations/directions for more accurate estimation.Such a strategy works well, but still incurs a large amount of additional computations.An alternative is to develop solutions that can estimate the grid bias directly, such as those proposed for off-grid DOA estimation in [24], [25], [26], [27], and [28].To the best of our knowledge, the off-grid source localization problem has not been addressed yet under the framework of sparsity maximization.In order to deal with the challenge, a two-step off-grid source localization method is also developed.In the first step, the source locations are roughly estimated with a coarser grid, whereas in the second step, their off-grid bias is estimated through an iterative process, which has a closed-form solution at each iteration.As a performance benchmark, the Cramér-Rao bound (CRB) is also derived.As demonstrated by computer simulations, compared with the grid refinement method, the proposed one can provide a better result with a lower computational complexity in terms of running time.
The main contributions can be summarized as follows.
1) The target localization problem with magnitude-only measurements is formulated into a group sparsitybased phase retrieval problem, where the measurements at all arrays are exploited simultaneously by enforcing the common spatial sparsity.2) To reduce the computational complexity, an offgrid estimation step is further introduced, where the off-grid biases of both x-axis and y-axis are approximated by Taylor series expansion with two variables.As a result, the proposed method does not require synchronization among different subarrays and it works on magnitude-only measurements, which leads to a robust solution where sensor phase errors have no effect on its performance, as also verified by simulations.

II. SIGNAL MODEL WITH DISTRIBUTED SENSOR ARRAYS
Consider K narrowband sources located at Cartesian coordinates L k (x k , y k ), k = 1, . .., K, impinging on D deployed sensor arrays with coordinates C d (x d , y d ), as shown is the steering matrix with its columns a(θ d,k ), k = 1, . .., K, being the corresponding steering vectors.When employing a UCA, a(θ d,k ) is given by [29] T (3) where γ m = 2π m/M d , λ is wavelength of the signals, r is radius of the circular array, and θ d,k denotes the arriving angle between the kth source and dth array, expressed as with y d,k = y k − y d , x d,k = x k − x d , and arctan2(•) being the inverse four-quadrant tangent operator.
For D distributed arrays, the measurements can be jointly expressed as where blkdiag{•} generates a block diagonal matrix from its entries, S = [s [1], . .., T .Note that magnitude-only measurements suffer from some ambiguities, and two of them affect the AOA estimation results: one is mirroring and the other is spatial shift [16], [17], [19], [21].The mirroring ambiguity refers to the phenomenon that signals arriving from −θ d,k will generate the measurements with the same magnitude, whereas for spatial shift ambiguity, it refers to that all estimated arriving angles at the array are phase shifted by a specific amount.However, as shown in [24], while applying a UCA structure, those ambiguities would not appear if the valid DOA range is limited within the range of 180 • ; therefore, UCAs are employed in this work.
By dividing the admissible area of interest into G x and G y grids along the x-axis and y-axis in the Cartesian coordinate system, separately, the overcomplete steering matrix of the dth array can be expressed as CORRESPONDENCE 9705 where is the angle between location (x g x , y g y ) and the dth array, with where y d,g x = y g x − y d and x d,g y = x g y − x d .Accordingly, the dth array measurements (1) are changed to where Sd = [s d [1], . .., sd [P]] is the sparse signal matrix with K rows corresponding to the targets being nonzero valued, and sd [p] represents the corresponding sparse signal vector at the pth snapshot with only K out of G entries being nonzero valued.Accordingly, a D d=1 M d × DG x G y steering matrix Ã covering all D sensor arrays can be constructed and the array measurements are given by where S = [s [1], . .., s[P]] is the joint signal matrix, with In practice, targets may not lay on a predefined grid.For the dth array, the steering vector corresponding to the real impinging angle θ d,k / ∈ {θ d,1 , . . ., θ d,G } can be approximated around its nearest on-grid angle θ d,g k by first-order Taylor expansion as [25], [26], [27] where β d,k denotes the off-grid bias and b(θ d,g k ) is the derivative of a(θ d,g k ).
For source localization, the true AOA at the dth array θ d,k is a function of two variables ( y d,g y , x d,g x ), determined by two unknown variables x k and y k , as presented in ( 6)- (8).As a result, while x k / ∈ {x 1 , . . ., x G x } and y k / ∈ {y 1 , . . ., y G y }, the steering vector of the dth array corresponding to the real source location can be approximated around its nearest on-grid location (x g k , y g k ) by first-order Taylor expansion of two variables as where β k,x and β k,y represent the off-grid biases with respect to x-axis and y-axis, separately, whereas b ∂y g are the partial derivatives of a(θ d,k ).Let x = [x 1 , . .., x K ] and y = [y 1 , . .., y K ] denote the true locations of K sources, and (8) can be approximated by the first-order approximation of two variables, given by and where Similarly, we have Bd,y = [b y (θ d,1 ), . .., b y (θ d,G )], y = diag(β y ), and where 2 , separately, where v x and v y are grid stepsizes with respect to x-axis and y-axis, respectively.
Noted that sources from an arbitrary grid point in the Cartesian coordinate system would share the same spatial support of Ãd , d = 1, . .., D, although the arriving angles with respect to different arrays are different.Therefore, the source localization problem can be formulated as a joint group sparsity-based optimization problem, given by min where with • 2,1 and • F represent l 2,1 norm and Frobenius norm, respectively, and ρ is the penalty term.The l 2,1 norm is defined as ŝg 2 (17) with ŝg being the gth row vector of Ŝ. Ŝ contains G = G x G y groups and the gth group, g ∈ {1, . .., G}, is a 1 × DP vector, consisting of gth row vectors of all Sd , d ∈ {1, .., D}.

B. Proposed Method
First, for the on-grid solution, i.e., the off-grid biases ˜ x and ˜ y are assumed to be zero, and the corresponding optimization problem in (15) The abovementioned formulation has the same form as those considered in group sparsity-based phase retrieval problem for DOA estimation, which can be solved by the algorithm ToyBar [21].
For the general off-grid case, instead of jointly estimate ˜ x , ˜ y , and S in (15), this problem is solved iteratively.In the first step, we employ (18) to find the rough grid locations of targets.In the second step, the PRIME technique [30] is employed, and the first term of (15) without off-grid bias at the pth snapshot is reformulated as where ãm represents the mth row of the steering matrix Ã, and z m [p] is the mth element of z where Re(•) represents the real part and the equality is met while s[p] = sq [p], it can be majorized as which can be formulated as where denotes the Hadamard product, sq is a known complex vector, and arg(•) represents the phase applied elementwise.Consider the off-grid biases, the PRIME technique is applied column by column to the first term of (15), and the corresponding objective function is changed to min with C = Z e jarg( ÃS e ) , where Se is estimated signals from step one, which is susceptible to global phase ambiguity.
After that, an iterative algorithm for estimating dictionary bias β is proposed.This method first estimates K nonzero rows of estimated signals Si where (•) † is the pseudoinverse operator, i ∈ {1, . .., I} is iteration index, and Āi K is the steering matrix with K columns corresponding to the estimated source locations, given by and θ i d,K is updated with x i and y i obtained from the previous iteration.x 0 and y 0 are initialized as K columns of Ãd , which corresponds to locations of estimated sources Se .As the biases βx and βy share the same support with Se , Bd,x and Bd,y are obtained, which are the submatrix of Bd,x and Bd,y .Similarly, the biases of corresponding source locations Once the PRIME technique is done, this process updates each off-grid bias βx and βy in (26) at a time.First, with βi y fixed to zero, βi x is pursued by solving Since Āi K and Bi x are block diagonal, ( 27) can be reformulated as where C d and Si d,K are the approximated measurements and estimated signals of the dth subarray.Dropping index i for simplicity, (28) can be approximated as [24], [25], [26] where tr(•) and Re(•) represent the trace and real part of its variable, separately.The optimal condition, βi x can be obtained by where (•) −1 is the inverse operator.While the off-grid bias with respect to x-axis is fixed, the off-grid bias of the y-axis β y can be obtained by (31) Similarly, we have (32) Note that, the noncoherent DOA estimation results still suffer from the global phase ambiguity, that is, where S r and A r represent the real signal and its corresponding real steering matrix, respectively, and φ is a global phase factor.When substituting (33) into ( 30) and ( 32), the global phase factor cancels, which means that the global phase ambiguity will not affect bias estimation in this step.With βi x and βi y , the steering matrix A d,K is updated as where β k,x and β k,y represent the bias of the kth source, and x 0 d,k and y 0 d,k are the initial locations obtained from the first step, i.e., corresponding locations of Se .Finally, the CORRESPONDENCE 9707 Algorithm 1: Algorithm Summary Input: Ã, Z, Initialization: Step 1: Estimate Se via existing group sparse phase retrieval algorithms Obtain Ā0 K , x 0 and y 0 Calculate C = Z e jarg( ÃS e ) .

Output localization results: x
estimated locations of the kth source can be obtained after I iterations.
The full algorithm is summarized in the Algorithm 1.

C. Cramér-Rao Bound
From ( 5), the probability density function is where the mth entry of Z at the pth snapshot and σ m is noise power.
Since the reconstructed signals are up to a global phase factor, for complex signals, the Fisher information matrix (FIM) would be singular [31].Thus, similar to [21], instead of estimating the phase information of signals, only phase differences between signals are considered.Assuming that the noise level at all arrays are identical, the unknown parameter vector of arriving angles, magnitude, phase difference, and noise level can be represented as is the phase of the kth signal of the dth sensor array at the pth snapshot.
For the deterministic model, the FIM F is defined as [32] The {i, j}th entry of F is given by [33] where −1 ( ) = 1 σ 2 m I M , I M is the identity matrix, (•) −1 is the matrix inverse operator, and μ( ) = |AS|.Since μ( ) is independent of the noise level, we have As the FIM is block diagonal, F σ m has no effect on the CRB result of locations.Thus, the location CRB can be determined by the inverse of F. Computing the derivatives of μ( ) with respect to , we have Denote |a m s| = (s H a H m a m s) (42) where (•) * is the complex conjugate operator, A m,d (k, :) is the kth row of A m,d , and A m,d (:, k) is the kth column of A m,d .The second block is expressed by Then, the third block is given by where stands for the Hadamard product.Then, F can be given by The CRB associated with locations of signals can be obtained by the inverse of F.

IV. SIMULATIONS
In this section, simulation results are provided to show the performance of the proposed methods in comparison with the existing sparsity-based on-grid method with coherent measurements (with both magnitude and phase information) [14].For on-grid methods, when grid refinement is employed, it is referred to as the refinement method, whereas for the off-grid method, the iteration number for the second step is 20.The sparse phase retrieval algorithm Toy-Bar in [21] is applied in the noncoherent scenario, and the number of iterations before stop is set to 300, with 30 random initializations used in order to find the global minimum.
The area of interest is set as [−20, 20] m along both x-axis and y-axis.In the initial step for both on-grid and off-grid methods, v x = v y = 2 m is used as the stepsize for constructing the overcomplete steering matrix Ã unless specified otherwise.In the refinement step, a new grid with stepsize 0. For the first set of simulations, the signal to noise ratio (SNR) is 20 dB.The spatial spectrum of the proposed noncoherent localization results is shown in Fig. 2, where Fig. 2(a) provides the result of the on-grid step, whereas Fig. 2(b) is for the extra off-grid step.It can be seen that the two sources have been identified successfully, but the off-grid step significantly increases the accuracy.
Next, performances of the three methods are evaluated with different SNR values ranging from 0 to 20 dB in terms of the root-mean-square error (RMSE) in the absence of phase error.The results are shown in Fig. 3, with each point obtained by averaging over 100 trials.It can be observed that, although all methods achieve more accurate results with increasing SNR, the method with full coherent measurements consistently outperforms those with magnitudeonly ones, especially when the noise level is high.This is not surprising since only magnitude information is used in the   noncoherent scenario, and the coherent method effectively provides a performance floor (in terms of RMSE, such as some kind of practical "CRB") when phase error is zero.Under the noncoherent scenario, the off-grid model clearly outperforms the on-grid model even when the on-grid model has a denser grid; moreover, the improvement achieved by a denser grid for the proposed off-grid method is not significant.In addition, the proposed noncoherent off-grid method with a 2-m stepsize has a better performance than the noncoherent grid-refinement method.
Then, we examine the performance in the presence of sensor phase errors.The array measurements with phase error is modeled as matrix with each entry being unit complex variable with a random phase term generated independently from the Gaussian distribution with standard derivation σ , representing the phase errors at the dth array.RMSE results are obtained with an average of 100 trials and the SNR is fixed at 20 dB.As shown in Fig. 4, the proposed noncoherent methods are not affected by phase errors, with a steady performance; on the other hand, the performance of the coherent method declines as the intensity of phase errors increases and becomes much worse than the noncoherent methods when the standard deviation σ of the phase error is larger than 0.15, and in this case, the proposed noncoherent methods are preferred.
To evaluate the performance for different number of snapshots, the RMSE results with phase error σ = 0 and σ = 0.2 are presented in Fig. 5, with SNR set to be 20 dB.Clearly, a larger number of snapshots yield more accurate results and the coherent method consistently outperforms CORRESPONDENCE 9709 the noncoherent one in the absence of phase errors.However, in the presence of phase errors, the coherent method suffers from a much larger RMSE.Finally, the computational complexity is compared in terms of running time with 100 snapshots, where the on-grid method took 226.11s, the off-grid one took 226.22s, and the refinement method took 395.84s,based on a computer with 4.2 GHz CPU i7-7700 K and 16 GB RAM.We can see that although the computational time of off-grid method is higher than the on-grid method, the time cost by the second step of the off-grid model is minimal, especially compared with the refinement method.

V. CONCLUSION
A novel source localization method with magnitudeonly measurements based on distributed array networks has been proposed.The noncoherent source localization problem was first formulated into a joint sparse phase retrieval form, and the l 2,1 norm is employed to enforce spatial sparsity, while to tackle the off-grid problem in the second step, dictionary bias is estimated through an iterative process with closed-form solutions at each iteration.As phase error at sensor arrays has no effect on the proposed solution, phase calibration is no longer required.Furthermore, the proposed off-grid method has provided more accurate results than the on-grid method with marginal computational cost, which is significantly less than the widely adopted grid-refinement approach.

Fig. 1 .
Fig. 1.Source localization geometry.in Fig. 1.The number of sensors at the dth sensor array is M d , and the corresponding noncoherent measurements are Z d = |A d S d | + N d (1) with S d = [s d [1], . .., s d [P]] and s d [p] = [s d,1 [p], . . ., s d,K [p]] T , where s d,k [p] represents the pth snapshot of the kth signal, N d is the M d × P random Gaussian noise vector at the dth array, | • |, and [•] T are the elementwise absolute value operator and matrix transpose operator, separately, and 2 m is formed around a distance of 1 m to either side of the estimated location from the initial step.There are D = 4 distributed sensor arrays placed at C 1 = (10, 40) m, C 2 = (30, 10) m,C 3 = (−80, 90) m, andC 4 = (−20, 40) m, whereas the off-grid locations for K = 2 sources are L 1 = (−10.5,−9.5) m and L 2 = (0.5, 12.5) m.The number of sensors at each distributed sensor array M d is set as 20, while the radius r of the UCAs is set as r = M d λ 2 2π , and P = 100 snapshots are collected unless specified, otherwise.