Investigation of Near-Field Source Localization Using Uniform Rectangular Array

: In ﬁfth-generation (5G) wireless communications, large-scale arrays pose a challenge to the accuracy of signal models based on the plane wavefront. In this paper, a novel method for 3D near-ﬁeld direction of arrival (DOA) estimation is proposed based on large-scale uniform rectangular array (URA). First, the near-ﬁeld signal model based on the vertical rectangular array and the delay phase shift of the received array is presented. Afterwards, the proposed method divides the complete parameters set into multiple-parameters subsets, and only estimates one of them in each iteration, leaving the others in the ﬁxed subset. As a result, we can obtain the maximum convergence rate of the deterministic maximum likelihood (DML) algorithm. Finally, the simulation results demonstrate that the root mean square errors (RMSEs) of the proposed algorithm are closer to the Cramer-Rao lower bounds and converge faster than those of the DML algorithm, conﬁrming its effectiveness and superiority.


Introduction
Source localization is a hot topic in numerous communication domains, such as wireless communications, radar, microphone arrays, and so on.In fifth-generation (5G) wireless communication systems, massive multiple-input multiple-output (MIMO) relies on the application of a large number of array elements at the base station, which is the key enabling technology [1].Massive MIMO can increase system throughput, spectrum efficiency, antiinterference capability, and high directivity.A large number of antennas not only improve the quality of communication but are also helpful in high-accuracy localization [1,2].Accurate DOA estimation is also crucial to establish the channel model of the massive MIMO system [1][2][3].On the basis of the plane-wave assumption, a number of high-resolution algorithms, such as the maximum likelihood estimation (MLE) method [4][5][6], the multiple signal classification (MUSIC) algorithm [7,8], the estimation of signal parameter via rotational invariance technique (ESPRIT) method [9][10][11], and some methods based on the optimization and extension of these algorithms [12][13][14][15][16][17][18][19][20], have been used for far-field source localization.
According to certain studies [21][22][23], the large antenna aperture and short wavelength allow signal sources to be within the Fresnel region (i.e., near-field) of the array aperture, thereby violating the plane wavefront assumption.In [22], the main focus was the LOS path azimuth angle variations for four different bands across the antenna arrays at different Tx locations.All of the findings show that the plane wavefront assumption is not valid, and the spherical wavefront should be considered (i.e., near-field).The capacity of a shortrange line of sight (LOS) MIMO system in 3D space was investigated in [24], and it was concluded that capacity statistics calculated using the spherical wave model (SWM) are more accurate and optimistic than those calculated using the plane wave model (PWM), with the difference becoming more obvious with the increase of the array aperture to transceiver distance ratio.As a result, the range information should be considered into the signal model [21,22].
For this purpose, many methods have been proposed to solve it, such as some improved methods [25,26] originating in far-field algorithms.In addition, a reduceddimension MUSIC algorithm based on uniform linear array (ULA) is proposed in [27] for near-field source localization, but this method is not suitable for a two-dimensional (2D) array.In [28], a novel near-field DOA estimation framework based on a deep complex-valued network is used for short-range massive MIMO, however the simulation results show that the proposed method has a better accuracy only in angle estimation.A search-free near-field source localization method is proposed in [29].However, the computational complexity of the parameter separation and polynomial rooting operations is high for multiple sources' localization.Some near-field localization algorithms based on the fourth-order cumulant are proposed in [30,31], but the implementation steps are found to be highly complex for a 2D array-based near-field model.
The maximum likelihood estimation (MLE) is an optimal estimator [5].The URA is one of the possible array configurations of massive MIMO and it can support more precise source location information by estimating azimuth, elevation, and range.However, lots of array elements make it hard for MLE to cope with the multi-dimensional optimum problem of the near-field model based on URA [28].The deterministic likelihood estimation algorithm based on the expectation-maximization (EM) method can solve this problem but suffers from the high computational cost.Therefore, a source localization method based on the near-field model is proposed in this paper, which utilizes second-order Taylor expansion to approximate the spherical wavefront [27,30].The parameter set of the DML algorithm is divided into multiple parameter subsets, and the maximum likelihood estimation (MLE) of the parameters is continuously iterated for each subset [32,33], which greatly reduces the computational complexity and has a faster parameter convergence rate.Numerical simulation results show that the proposed method can perform better than the DML algorithm in adverse cases, such as low signal-to-noise ratio (SNR) or correlated signals.
The paper is structured as follows.In Section 2, the near-field model based on the vertical rectangular array and the delay phase shift of the received array is presented.Section 3 introduces the DML algorithm principle.In Section 4, the proposed method is provided in detail.Numerical simulation results show the superiority of the proposed method in Section 5. Finally, the summary is concluded in Section 6.

Near-Field Signal Model
Consider a near-field scenario at the receiver, where K narrowband source signals are impinging on the uniform rectangular array located in the XOZ plane with (2M + 1)(2N + 1) array elements and with inter-element spacing, d, along each axis, as depicted in Figure 1.In this case, its central array element can be considered as the phase offset reference point and the coordinate origin.The received signals at ) array element can be expressed as [28,34,35] where  () is the -th ( = 1,2, ⋯ , ) signal received at the reference point,  ,  , and  are the azimuth, the elevation, and the range, respectively, for the -th signal,  presents the wavelength of the signals, , , denotes the -th source's amplitude ratio between the (m,n)-th sensor and the reference point,  , , is the distance between the (m,n)-th sensor and the -th source,  is the distance between the reference element and the -th source, and  , () is the additive white Gaussian noise (AWGN) at the (m,n)-th sensor element.
. In the near-field model, when the -th signal impinges on the array,  , , can be calculated as: Where Δs is the wave path difference between the (m,n)-th sensor and the reference point in the far-field.According to the geometry, Δs =  cos  sin  +  cos  .The wavefront's curvature is no more planar when sources are located close to the array, i.e., in the near-field or Fresnel region ( < 2 /) [29], where  is the array aperture (the long diagonal length of the URA).If we retain terms up to the second power of , then its multiplier after using the binomial expansion theorem can be written as: Here, , , = 2  ( cos  sin  +  cos  ) (8) Additionally, where s k (t) is the k-th (k = 1, 2, • • • , K) signal received at the reference point, θ k , ϕ k , and r k are the azimuth, the elevation, and the range, respectively, for the k-th signal, λ presents the wavelength of the signals, r k r m,n,k denotes the k-th source's amplitude ratio between the (m,n)-th sensor and the reference point, r m,n,k is the distance between the (m,n)-th sensor and the k-th source, r k is the distance between the reference element and the k-th source, and z m,n (t) is the additive white Gaussian noise (AWGN) at the (m,n)-th sensor element.
In the near-field model, when the k-th signal impinges on the array, r m,n,k can be calculated as: Where ∆s is the wave path difference between the (m,n)-th sensor and the reference point in the far-field.According to the geometry, ∆s = md cos ϕ k sin θ k + nd cos θ k .
The wavefront's curvature is no more planar when sources are located close to the array, i.e., in the near-field or Fresnel region (r F < 2D 2 /λ) [29], where D is the array aperture (the long diagonal length of the URA).If we retain terms up to the second power of d r k , then its multiplier after using the binomial expansion theorem can be written as: Here, Note that the error caused by the Fresnel approximation ( 4) is significant if r F ≤ 0.62 D 3 /λ 1/2 .In matrix form, (1) can be expressed as: where, ×K is the array manifold matrix, and can be written as: with

and [•]
T de- notes the transpose.
In this paper, the following assumptions must be held: (a) Z(t) is the zero-mean, temporally, and spatially additive white Gaussian noise with covariance matrix σ 2 I, which is independent of the incident signals, where I is the identity matrix.(b) The source signals are unknown but deterministic signals.
When the number of snapshots is L, the estimated problem of this paper is to determine the DOAs {θ and S from these snapshots, §, where:

Deterministic Maximum Likelihood Estimation
According to the MLE method, the unknown parameter set {θ, ϕ, r} and S can be solved by the following optimization expression [25]: The deterministic maximum likelihood estimation can go a step further to solve ( 14), as follows [25]: First Step: For the given {θ, ϕ, r}, the estimated values of Ŝ(t) that minimize ( 14) can be calculated as: Second Step: Replace Ŝ(t) with S(t) in ( 14) and obtain the estimation of {θ, ϕ, r}: where  [25].However, in the second step, the maximization of ( 16) with respect to {θ, ϕ, r} is a complex multi-parameter optimization problem.Therefore, the EM method can be used to solve such problem [6].For the EM solution to our problem, both complete data and incomplete data are utilized.In the proposed model, the received superimposed signals are disturbed by additive white Gaussian noise, so the mixed signals which are composed of signal and noise can be observed at the receiver as the incomplete data [25].
It is almost impossible for us to estimate these parameters based on the incomplete data, Y(t), but the log-likelihood function of the complete data is simple to maximize (16) and the complete data log-likelihood function can be estimated from the incomplete data without any difficulty.Therefore, we can explain {θ k , ϕ k , r k } K k=1 by constructing natural complete data that are affected by partial noise.The design process can be summarized as: where Gaussian noises with the zero-mean whose variance is The relationship between the incomplete data, Y(t), and the complete data, y k (t), is expressed as: Define: As a result, the log-likelihood function of the complete data, †, is calculated as [6]: The EM method consists of the two-step iterative procedure, i.e., E-step (expectation) and M-step (maximization), and it utilizes (22) to complete the estimation of unknown source parameters.When we are at the (p + 1)-th iteration, the two-step process is as follows [6, 25,32,33]: 1. E-step: Calculate the conditional expectation of the complete data log-likelihood according to the previous estimated results of {θ p , ϕ p , r p } and S p :

M-step: Maximize
. By the simplifications of ( 23), the conditional mean can be obtained as [25]: To the right of (24), the first term a θ P k , ϕ p k , r P k s p k (t) is the known signal estimated at the p − th iteration, while the second term of ( 24) is equivalent to β k times of the overall noise vector.According to the EM method, is orthogonal to the signal subspace which depends on A(θ p , ϕ p , r p ) [25].
In the M-step, after inserting y p+1 k (t) into (22), the maximization of the complete data log-likelihood regarding {θ k , ϕ k , r k } and s k (t) can be given as: where Kp H is the covariance of † k .In addition, ref. [32] shows that the convergence rate of the EM method is inversely proportional to the Fisher information of the constructed complete data.If all unknown parameters are asked to update at the same time in one iteration, only each component of the complete data contains the same noise power, i.e., K .The complete data can obtain the fastest convergence rate because they contain the least Fisher information.Notably, the EM method does not specify all unknown parameters to be updated at the same time in one iteration.

Parameter Sequential Updating Strategy Based on SAGE
As a result, the Space-Alternating Generalized Expectation-Maximization (SAGE) method can be used to improve the parameter sequential updating of the EM method, i.e., when we estimate the unknown parameters of the k-th signal, assuming that all parameters of other signals are known and fixed [32,33], we only construct the complete data for the parameters of the k-th signal, which reduces the Fisher information quantity of the complete data.
Afterwards, we divide {y k (t)} K k=1 into I groups, and each group is composed of K i signals (i = 1, • • • , I).The signals' parameters in the other groups are assumed as known when the parameters of the i-th group are estimated.The selected i-th group is "complete data", and the dimensions of the complete data are reduced to K i when the noise power of each component of the i-th group is equal to 1 K i times to the noise power of the observation data, †, and it has the lowest Fisher information.
The Fisher information of the complete data after grouping is lower than the Fisher information of the complete data before grouping due to K i < K.As a result of the higher step size of parameter updating, it can reach a stable point faster.When I = 1, the SAGE method is equivalent to the EM method, and if I = K, the iteration process can achieve the fastest convergence rate.Considering all of that, {β k (t)} K k=1 can be defined as follows [32]: We further divide {θ, ϕ, r} into K subsets, where each subset is: To reduce the computing complexity, we replace the 3D optimization approach for calculating the MLE of one source signal's parameters with three one-dimensional (1D) methods, where every parameter is estimated separately [32].It is attained by splitting each subset {θ k , ϕ k , r k } into three subsets, i.e., {θ k }, {ϕ k }, and {r k }.For each of these subsets, ( 17) is still applicable.For re-estimating the parameter in the three subsets, the following updating approach can be utilized to estimate the series connection of the three SAGE rounds: The initialization is very important to the global convergence of the proposed algorithm.Hence, we present an effective initialization method after we define Θ k = {θ k , ϕ k , r k }, and the parameter vector of the SAGE approach is: Firstly, we can regard Y(t) as a single source before starting the initialization steps.We therefore obtain: Next, we can obtain s 0 1 by ( 15), then: We repeat these steps until we obtain the initial values of all signals.Finally, Considering all the above, after the i-th iteration, we assume that the estimation of all parameters set Θ is

10: end
This shows that the computational cost of one iteration step of the DML algorithm is identical to that of one iteration cycle of the proposed algorithm without ( 29)- (31).In addition, the SAGE-based parameter updating strategy can improve the convergence rate, resulting in less iteration.

Simulation Results
The validity of the proposed method is tested using the RMSE (root mean squared error) of 600 Monte Carlo simulations under different cases: where r k is the real value of the range parameter, and rk,l is the estimated value of r k in the l-th test.The RMSEs of other parameters, such as θ k and ϕ k , have a similar form as that of (36).Note that the estimated range is the rk,l -to-wavelength ratio.
In the proposed model, the CRLB is given by the following expression [36,37]: with: where {•} denotes the real part of the given argument, denotes the Schur-Hadamard product, and denotes the Kronecker product.J P×P is the P × P matrix of one, and P is the number of the unknown parameters.

CASE 1
Consider a URA with M = N = 4, and the interval between two elements is set at d = λ/2.Therefore, the near-field region is r F ∈ [8.35λ, 64λ].Two uncorrelated sources are located at 20 • , 100 • , 9λ and 120 • , 70 • , 12λ , respectively.Furthermore, the number of snapshots is 50, and SNR varies from −20 to 30 dB.The RMSEs of the proposed method are smaller than those of the DML algorithm at different SNRs, as shown in Figures 2-4.When the sources are well-separated and the signals are uncorrelated, they are both extremely close to the theoretical bound for high SNRs.Moreover, in comparison to the DML, the RMSEs of the proposed method perform better at a low SNR.where  is the real value of the range parameter, and ̂ , is the estimated value of  in the l-th test.The RMSEs of other parameters, such as  and  , have a similar form as that of (36).Note that the estimated range is the ̂ , -to-wavelength ratio.
In the proposed model, the CRLB is given by the following expression [36,37]: with: where ℜ • denotes the real part of the given argument, ⊙ denotes the Schur-Hadamard product, and ⨂ denotes the Kronecker product. is the   matrix of one, and  is the number of the unknown parameters.

CASE 1
Consider a URA with   4, and the interval between two elements is set at  /2.Therefore, the near-field region is  ∈ 8.35, 64 .Two uncorrelated sources are located at 20 °, 100 °, 9 and 120 °, 70 °, 12 , respectively.Furthermore, the number of snapshots is 50, and SNR varies from −20 to 30 dB.The RMSEs of the proposed method are smaller than those of the DML algorithm at different SNRs, as shown in Figures 2-4.When the sources are well-separated and the signals are uncorrelated, they are both extremely close to the theoretical bound for high SNRs.Moreover, in comparison to the DML, the RMSEs of the proposed method perform better at a low SNR.

CASE 2
Consider a URA with   4, and the interval between two sensors is set at  /2 .Therefore, the near-field region is  ∈ 8.35, 64 .Two correlated narrowband sources are located at 30 °, 40 °, 8.7 and 60 °, 50 °, 12.5 , respectively.The correlation coefficient between signals is set as 0.8762.The number of snapshots is  50, and the SNR is −5 dB.As shown in Figure 5, in adverse environments including low SNR and high correlation, because of the faster convergence rate, the RMSEs of the proposed method converge in 4-10 iterations, but those of the DML algorithm converge in 8-18 iterations.Besides, we also conclude that the parameter estimation error of the proposed method is smaller than that of the DML in general.

CASE 2
Consider a URA with   4, and the interval between two sensors is set at  /2 .Therefore, the near-field region is  ∈ 8.35, 64 .Two correlated narrowband sources are located at 30 °, 40 °, 8.7 and 60 °, 50 °, 12.5 , respectively.The correlation coefficient between signals is set as 0.8762.The number of snapshots is  50, and the SNR is −5 dB.As shown in Figure 5, in adverse environments including low SNR and high correlation, because of the faster convergence rate, the RMSEs of the proposed method converge in 4-10 iterations, but those of the DML algorithm converge in 8-18 iterations.Besides, we also conclude that the parameter estimation error of the proposed method is smaller than that of the DML in general.

CASE 2
Consider a URA with M = N = 4, and the interval between two sensors is set at d = λ/2.Therefore, the near-field region is r F ∈ [8.35λ, 64λ].Two correlated narrowband sources are located at 30 • , 40 • , 8.7λ and 60 • , 50 • , 12.5λ , respectively.The correlation coefficient between signals is set as 0.8762.The number of snapshots is L = 50, and the SNR is −5 Db.As shown in Figure 5, in adverse environments including low SNR and high correlation, because of the faster convergence rate, the RMSEs of the proposed method converge in 4-10 iterations, but those of the DML algorithm converge in 8-18 iterations.Besides, we also conclude that the parameter estimation error of the proposed method is smaller than that of the DML in general.

CASE 3
Consider a URA with M = N = 4, and the interval between two elements is set as d = λ/2.Therefore, the near-field region is r F ∈ [8.35λ, 64λ].In the first group, two uncorrelated sources are located at 30 • , 40 • , 9λ and 30 • , 45 • , 10λ , respectively.In the second group, two uncorrelated sources are located at 30 • , 40 • , 9λ and 30 • , 100 • , 10λ , respectively.Moreover, the number of snapshots is L = 100, the SNR is 20 dB, and the iteration varies from 1 to 30. Figure 6 illustrates that the convergence rate is faster if the sources are well-separated, and the RMSEs of the well-separated sources are smaller than those of closely spaced sources.Figure 6b shows that the convergence rate will be affected if the sources are located closer to each other.
Consider a URA with  =  = 4, and the interval between two elements is set as  = /2.Therefore, the near-field region is  ∈ [8.35, 64].In the first group, two uncorrelated sources are located at (30 °, 40 °, 9) and (30 °, 45 °, 10), respectively.In the second group, two uncorrelated sources are located at (30 °, 40 °, 9) and (30 °, 100 °, 10), respectively.Moreover, the number of snapshots is  = 100, the SNR is 20 dB, and the iteration varies from 1 to 30. Figure 6 illustrates that the convergence rate is faster if the sources are well-separated, and the RMSEs of the well-separated sources are smaller than those of closely spaced sources.Figure 6b shows that the convergence rate will be affected if the sources are located closer to each other.

CASE 4
Consider a URA with M = N = 4, and the interval between two elements is set as d = λ/2.Therefore, the near-field region is r F ∈ [8.35λ, 64λ], and two uncorrelated sources are located at 20 • , 120 • , 9λ and 30 • , 40 • , 10λ , respectively.Moreover, the number of snapshots is L = 100, the SNR is 10 dB, and the iteration varies from 1 to 30.In addition, consider a URA with M = N = 4, and the interval between two elements is set as d = λ/2.Therefore, the near-field region is r F ∈ [8.4λ, 64λ], and two coherent sources are located at 20 • , 120 • , 9λ and 30 • , 40 • , 10λ , respectively.The number of snapshots is L = 100, the SNR is 10 dB, and the iteration varies from 1 to 30.As shown in Figure 7, the RMSEs of all parameters can converge in 8-10 iterations whether they have uncorrelated signals or coherent signals, and the proposed method is not very sensitive to source correlation.
Consider a URA with  =  = 4, and the interval between two elements is set as  = /2.Therefore, the near-field region is  ∈ [8.35, 64], and two uncorrelated sources are located at (20 °, 120 °, 9) and (30 °, 40 °, 10), respectively.Moreover, the number of snapshots is  = 100, the SNR is 10 dB, and the iteration varies from 1 to 30.In addition, consider a URA with  =  = 4, and the interval between two elements is set as  = /2.Therefore, the near-field region is  ∈ [8.4, 64], and two coherent sources are located at (20 °, 120 °, 9) and (30 °, 40 °, 10), respectively.The number of snapshots is  = 100, the SNR is 10 dB, and the iteration varies from 1 to 30.As shown in Figure 7, the RMSEs of all parameters can converge in 8-10 iterations whether they have uncorrelated signals or coherent signals, and the proposed method is not very sensitive to source correlation.

Computational Complexity
As described in Table 1, we only focus on the computational complexity of the DML and the proposed method in (29)- (31).Here, ∆θ, ∆ϕ, and ∆r represent the search steps of elevation, azimuth, and range, respectively.L is the number of snapshots, and (2M + 1)(2N + 1) is the number of array elements.By looking at the complexity histogram as shown in Figure 8, we can gain a better sense of the difference in computing cost between the two algorithms.and the proposed method in ( 29)- (31).Here, ∆, ∆, and Δr represent the search steps of elevation, azimuth, and range, respectively. is the number of snapshots, and 2 1 2 1 is the number of array elements.By looking at the complexity histogram as shown in Figure 8, we can gain a better sense of the difference in computing cost between the two algorithms.The run-time is presented in Table 2 and the CPU is AMD 5800U.It can be seen that the search step has a significant impact on the computational cost, and the proposed algorithm has a smaller computational cost than the DML algorithm.The run-time is presented in Table 2 and the CPU is AMD 5800U.It can be seen that the search step has a significant impact on the computational cost, and the proposed algorithm has a smaller computational cost than the DML algorithm.

Conclusions
In this paper, we proposed a computationally efficient near-field source localization method based on URA.The proposed method does not involve 3D optimization, polynomial rooting operations, or high-order cumulants' calculations.The Monte Carlo simulations demonstrated that the performance of the proposed method compared to the DML algorithm was closer to the CRLBs at a high SNR and performed better in adverse cases such as low SNR or correlated signals.Moreover, we observed that the convergence rate was affected if the sources were located closer to each other, and we also conclude that the proposed method had a smaller computational cost than the DML algorithm.Future work will concentrate on the moving and mixed near-to far-field sources.

Figure 1 .
Figure 1.Spherical wavefront propagation model in the near-field condition.

Figure 1 .
Figure 1.Spherical wavefront propagation model in the near-field condition.

Figure 2 .
Figure 2. RMSE and CRLB versus SNR of the estimated parameter .

Figure 2 .
Figure 2. RMSE and CRLB versus SNR of the estimated parameter θ.

Figure 3 .
Figure 3. RMSE and CRLB versus SNR of the estimated parameter .

Figure 4 .
Figure 4. RMSE and CRLB versus SNR of the estimated parameter .

Figure 4 .
Figure 4. RMSE and CRLB versus SNR of the estimated parameter .

Figure 4 .
Figure 4. RMSE and CRLB versus SNR of the estimated parameter r.
(a) The RMSE of elevation versus iterations in adverse conditions.(b) The RMSE of azimuth versus iterations in adverse conditions.(c) The RMSE of range versus iterations in adverse conditions.

Figure 5 .
Figure 5. Convergence comparison of the two algorithms.Figure 5. Convergence comparison of the two algorithms.

Figure 5 .
Figure 5. Convergence comparison of the two algorithms.Figure 5. Convergence comparison of the two algorithms.
(a) The RMSE of elevation versus iterations in different source spacing.(b) The RMSE of azimuth versus iterations in different source spacing.(c) The RMSE of range versus iterations in different source spacing.

Figure 6 .
Figure 6.Convergence comparison of the source spacing.Figure 6. Convergence comparison of the source spacing.

Figure 6 .
Figure 6.Convergence comparison of the source spacing.Figure 6. Convergence comparison of the source spacing.
(a) The RMSE of azimuth versus iterations in different correlations.(b) The RMSE of elevation versus iterations in different correlations.(c) The RMSE of range versus iterations in different correlations.

Figure 7 .
Figure 7. Convergence comparison of coherent signals.Figure 7. Convergence comparison of coherent signals.

Figure 7 .
Figure 7. Convergence comparison of coherent signals.Figure 7. Convergence comparison of coherent signals.

Figure 8 .
Figure 8. Complexity comparison of two algorithms.

Figure 8 .
Figure 8. Complexity comparison of two algorithms.

Table 1 .
Complexity comparison of the two algorithms.

Table 2 .
Run-time of the two algorithms.

Table 2 .
Run-time of the two algorithms.